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SONIC BOOM PROPAGATION IN A STRATIFIED ATMOSPHERE, 
WITH COMPUTER PROGRAM 


Wallace D. Hayes*, Rudolph C . mef< eli, and H E. 
Aeronautical Research Associates of Princeton, 


Kulsrud 
Inc . 


SUMMARY 


An analysis is presented of the propagation of sonic boom m 

ysrsss ."K^srs^J* 

Reauired inputs include atmospheric properties and horizontal 
winds Is functions of altitude, information on the flight path of 

?aSa£?e! The output includes midfield pressure signatures at any 
altitude . 

Results from sample calculations are presented and discussed. 


INTRODUCTION 


Sonic booms have become of prime importance in the design and 

Far Her algorithms for sonic boom have used various unjusti 
fied simplifying°assumptions . A basic aim of Present algo- 
rithm has been to avoid these assumptions as faras po ®^ b ^® e ^ t 
to extend the cases which could be considered. Thu , P 

algorithm includes the following features: 


(1) The inclusion of maneuvering aircraft in a sonic boom 
pressure calculation; 

(2) An appropriate ray-tube area calculation based on linear 
geometric acoustics; 

♦Professor of "Aerospace Sciences, Princeton University 



^ the ° f cora P let e signatures , without 

”^o" leld • a ^? Um ? t:LOnS;, obta ined through the use of an 
age variable m the calculation of nonlinear effects. 

^h P ^ 6S ? nt Algorithm assumes a horizontally stratified atmosohere 
with horizontal winds but without turbulence This limiting 
assumption corresponds to the case of greatest practical interest 
and considerably simplifies the calculation. praci;ical ^tereso 

The analysis is largely a rational synthesis of existing 

development SCr SDecific literature ’ with some new theoretical 
ueveiopment . .Specific references are cited in the bodv of iho 

laiM on * r prijac d pal new theoretical development is in the calcu- 
lation of ray-tube area. The analysis is also new in the care fin 
piecing together of a number of calculations, principally in the 
r ?if'^ lon oP bb e wave system and rays issuing from the aircraff 
stratified e at ySte h ^ r ?? S properl y describing propagation in the 

Sonsideiatio^ o? P a 6 Si^ w j nds ' Thls relation requires the 

consiaeration of a galilean transformation connecting* a innai 

dinate system with the fixed coordinate systeS? 8 ° Cal COOr ' 
A note of caution at this point may be in order Althonp-h 

?Seor?Ss yS it S desc ^ bad aa largely a syntSsS'of eSsWng 

heories, it should be pointed out that not all these theories mav 
be familiar to all workers in the field of sonic boom In order J 
to make the analysis feasible, the concept of galileaA invariance 
has been brought in from the subject of mechanics and a number nf 

pr opaea t i on Ve in frJm the generSf ?heo^o? wave^ ° f 

soS?cfJ A*numSe^ P n? i? en llterature ls diffused through many 
erent from ihaff? f b sac papers were written in contexts diff- 

ouTLS^i^ a syn?iesS°oi- "f 

sonic^boom^ 

tv ( |fS ibal computer Program has been written in ASA FORTRAN 
f SOme llteral text enclosed in asterisks) with flexibi- 

varLL^f 11 a ^ m * The Program is designed to be usable on a wide 
variety of modern computers and to be applicable to a variof-U^? 

problems It was developed using an IBM-H30, Model 2B and^hen 
modified for and operated with a CDC-66OO. The program mav be 
ered to accommodate the operating system constraints of a 

dtslmatiL? 0m ir™v fe a?«o S S SlmPle Changes ln input-output unit 

program structure from subprogram 3 linage tr^ai^rragraiTlirkLr 11 

rh^ e - S R Ci £ 10atlons of ln P ut “<3 units „ h t ° s 

phere is to be specified, and in how certain maneuver time deX 
tives are to be obtained from Input data. "‘ euve r time denva- 
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This report is accordingly divided into developmental the 

SS;fo^;Se S o^r S -riding and 

t ? ? rSafoS^far^ 1 " 

tion of the theory , with accent o P t he course of the analysis, 

motivation underlying the anaiys . llcatlons q n the computer 

rr 1 o e S ra^! ate Brs?^s a cur^h 1 t U ?Se?enc^, a ?S^ a are s« historical 
notes appended. 

The second part, consisting of COMPUTER PROGRAM “d OOMPUTA- 

oS^pSt^istSls’are included! aid typical computation results are 
presented . 

SYMBOLS 


This section includes symbols i ^ e ^® e “ al |§RTRAN X !!!bo? 

Shira^eSp^S the^compuLr^program^are identified in 

Table 1. 


Symbol 

a 

A 


A(x q ) 


n 


o 


"P 


'D 


L 


page No 


speed of sound (eq. 6) 

21 

ray-tube area cut by horizontal plane 

34 

(eqs. 19 > 26 ) 

area distribution of slender body 

43 

group or ray velocity (eq. 12) 

29 

normal phase velocity (eq. 11) 

28 

Snell's law invariant (eqs. 14, l6) 

29 

specific heat at constant pressure 

65 

wave drag coefficient (eq. 36 ) 

49 

drag coefficient per unit azimuth angle 

(eq. 36 ) 49 

lift coefficient (eq. 4) 

18 
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18 


°T C D 
D 


f ,f 
Y z 

F 

F. 


F 


f 


G 


' 1 , 2,3 


H 
I. 

k 

K R 

i ,m,n 

£(x 1 ) 

L 

L 

L a 
M 

n 


n Ik n T 

N(t) 


P 

q 

o 

r 

r 


net axial force coefficient (eq. 3) 


drag (eq. 3) 

line force distributions 


F function for aircraft signatures (eq. 33) 
input F-function 

F-function conversion factor (eq. 34) 
gravitational acceleration (eq. l) 
altitude of ground above sea level 
integrals used in calculation of A (eq. 24 ) 
heat conduction coefficient (eq. 58) 
reflection factor 

direction cosines of initial wave normal (eq.8) 
equivalent line force distribution 
lift (eq. 4) 

distance along aircraft axis (local phase) 
aircraft reference length 


46 

46 

48 

16 

16 

39 

63 

64 
25 
45 
18 
48 
48 


Mach number of aircraft, v/a 

unit vector normal to wave front 

normal and axial load factors (eqs. 3, 4) 

reduced longitudinal kinematic viscosity 
(eq. 60 ) J 


21 

28 

18 

66 


pressure (eqs. 1, 39) 
perturbation velocity (eq. 39) 
dynamic pressure (eq. 3) 

cylindrical radius in local coordinates 
position vector (eq. 12) 


16 

51 

18 

44 

29 
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r 1 

horizontal position vector 

36 

R 

gas constant; hyperbolic radius 

44 

s 

distance normal to wave front (local phase) 

56 

Ref 

reference wing area for force coefficients 
( eqs . 3, *0 

18 

S(x 1 , 0 r ) 

area distribution of equivalent body of 
revolution 

45 

s o U) 

integral of Vg(|) (eq. 54) 

63 

1 — l 

integral of Vg(^px) (eq. 51 ) 

58 

t(x 0 ,y 0 ) 

t 

wing thickness 

time along ray (eq. l 8 ) 

44 

26 

t a 

time along aircraft trajectory (eq. 2 ) 

18 

T 

absolute temperature 

l 6 

T 

thrust (eq. 3 ) 

18 

u 

wind velocity ( _u x ^ _u y^°) ( eq * 2 ) 

16 

U ,11, 

n t 

minus components of u in (x-, ,y n ) coordinates 
(eq. 15) 

31 

V 

aircraft velocity relative to atmosphere, Ma 
(eq. 2 ) 

17 

V E 

measure of signal invariant on kinematic ray 

(eqs. 40, 45) 

52 

W 

weight of aircraft (eqs. 3> 4) 

18 

i'x,y,z 

fixed coordinate system; east, north, and 
above ground, respectively 

16 

<j x’ ,y ' , -z 

t 

coordinate system aligned with aircraft 
velocity 

22 

[x 1 ,y 1 , -z 

coordinate system aligned with wave normal 
(eq. 19) 

31 
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x,y,z 


x ,y ,z 
o ’ o o 


X-, 


P 

P 

7 

A, 6 

4 

9 


M- 

v 


I 


1 


P 

T 


0 

0 

0 


a 

r 


$ 

f 


local coordinates near aircraft 43 

dummy coordinates near aircraft 43 

axial coordinate for equivalent body of 45 

revolution 

Prandtl-Glauert parameter, (M 2 - l) 1//2 43 

inverse of atmospheric scale height 59 

aircraft climb angle (eq. 2) 17 

ratio of specific heats 16 

perturbation from undisturbed value 39,28 

wind heading angle (whence wind blows) 16 

inclination angle of n below horizontal 25 

( eqs . 10, 17) 

Mach angle, sin -1 ( 1 /M) = tan -1 ( 1 /P) 22 

shear and dilatational viscosities (eq. 58) 65 

heading angle of wave normal n (eq. 9) 25 

linear phase variable (time) (eq. 4l) 51 

actual phase variable (time) (eq. 49) 56 

atmospheric density (eq. l) 16 

age (eq. 46) 57 

azimuth angle of wave normal from vertical plane 22 

aircraft bank angle 18 

azimuth angle of wave normal relative to 20 

aircraft 

local perturbation velocity potential (eq.30) 43 

heading angle of aircraft (eq. 2) 17 
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Subscripts 

a aircraft 

o initial value at time of emission of a ray 

from aircraft 

w wing 


17 

32 

45 


Vector Components 


Vector 


Name 


Symbol 


Position 


Coordinate Systems 


(x,y,z) 


x,y,z 




' ,y' , -2 


( x i 5y i ?-z ) 


x i' y i' ~ z 


Horizontal unit. 


1,0,0 

east 

i 

Horizontal unit. 


0,1,0 

north 

7 

Vertical unit 

k 

0,0,1 

Horizontal pos- 


x,y,0 

ition 

r 1 

Aircraft 

V 


velocity 

i 



Initial wave 



normal 

n o 


Wave normal 

n 


j 

IWind 

u 

1 

£ 

><1 

1 

£ 

<<! 

O 

Horizontal unit, 


sin f ,cos 

propagation 

n ! 

Inverse phase 

-1- 


velocity 

c n 
n 


Ray or group 



velocity 

c 



sin f , -cos ^ , 0 

sin v , -cos v ,0 ■ 

cos if, sin f, 0 | 

cos v, sinv ,0 

1 — l 

1 

o 

•\ 

o 

o 

0 

1 

H 

t l ~ 

X ,y , o 

x 1 ,y 1 ,o 

V cos 7 , 0 , 


-V sin 7 


f ,m,n 

cos 0 . 0 , 
sin 0 ° 


o 


I 

| c os 0 , 0 , sin 0 

i 


if, o 


1 , 0,0 

c^O^'Han 0 


a cos 6 - u 
-u^, a sin 
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THEORETICAL ANALYSIS 


General Description 

The intent of this section is to present a guide to the 
theoretical analysis which will be developed in this chapter. 

This guide is presented in several subsections. The first gives 
a brief description of the nature of sonic boom theory. ^The next 
three discuss certain basic concepts of geometric acoustics, with 
one purpose being that of explaining the basis for the assumption 
of steady ray geometry on which the entire analysis is based. The 
last describes the detailed analysis In digest form, essentially 
section by section. 

The nature of sonic boom theory .- Sonic boom is an acoustic 
phenomenon The appropriate theory for sonic boom propagation 
In acoustic theory with both simplifications and complicatrons 
which do not normally appear in acoustic theory With cert 

IfoSttir’anSoglSn^feLel^LILS: 5he ?heor? of geometric 

acoustics' is valid in an asymptotic sense when the wave length 1 

small compared with characteristic macroscopic scales of th 

problem. Such macroscopic scales include the of the atmos- 

of the wave fronts and the scale height (e.g p/p 6 ) of the aumo 
nhere (Symbols are defined in the list of SYMBOLS ) Geometric 
ac oust ic s is invalid in the region . i near _ the 

separate treatment is needed to obtain initial conditions 
propagated signal. 

Standard acoustic theories are linear. In 

nonlinear effects are locally very weak, but they have a 
nnnripeli ^ible cumulative effect during propagation over large 
d?stencel ^ ?he cumSiative nonlinear effect comprises distortion 

fesSibffSnL*^^ acoustics, 

theoretical approach to sonic boom which is here developed 
detail may be found in reference 1. 

Tn anv sonic boom theory with the generality of _ the theory 
presented here, the concepts of galilean. transformations an o 
phase necessarily appear. These concepts are discussed bel w. 

Coordinate systems and galilea n transformations .- Two ^prin 
cipal coordinate systems are required in “V theory One is n 
unaccelerated coordinate system fixed relative to the ground 

In ^Merited coordinate system aligned with the aircraft 
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of f?5 a *If and mo Y ang th aircraft velocity at the instant 
interest. .The flow near the aircraft is conveniently described 
m this coordinate system. These two coordinate systems are 
related through a galilean transformation. A galilean transform- 
er 1011 is a transformation from one unaccelerated coordinate system 
mov t?f rel ative to the first at a constant velocity. A 
quantity is galilean invariant if it does not change under a 
galilean transformation. 

ane particular step of the analysis the consideration of 
the galilean _ transformation is inescapable. This step appears 
en the variables describing the (local) flow near the aircraft 
are transformed . into the appropriate variables describing the 
(glo al) acoustic propagation in the coordinate system fixed rela- 
ll ve It 6 | round - In. this critical step we shall avoid going 
r ?^ formad details of the galilean transformation. Instead 

we identify corresponding variables which are inherently galilean 

of V the a nrnbiJ relatlng t J ese to both td e local and global variables 
of the problem, we are able to connect the local with the global 

variables This stratagem simplifies this critical step consid- 
erably and, . m effect, _ accomplishes the inescapable galilean 

ln a rel atively easy way. No other feasible way of 
relating the local and global variables was discovered. 

^ In + t ?^- s r ep°rt we are concerned primarily with the case of a 
nt ? 1:L I st ^ a ^ afied atmosphere with winds. Such an atmosphere 
remains horizontally . stratified under a horizontal galilean trans- 
formation, one m which the relative velocity is horizontal. 

Hence, any theory for this case must be invariant under such a 
transformation. This property has been used in the development of 
the analysis presented here to check it for algebraic consistency. 

. S e n.eral, consideration of which variables are galilean in- 
variant was of great help in the development of the analysis 

J hi I f° port ‘. A quantity which is galilean invariant 
If ande P enden 't °P Ihe choice of coordinate system. It is found 
that trie analysis is simpler, both algebraically and conceptually, 
when such variables are chosen to describe the solution. Thus th4 
consideration. of galilean invariance has guided the general course 
o the analysis and the specific choice of variables used. 

In this report, we mention in a number of places whether 
particular variables are or are not galilean invariant. Except in 
ie cr itical step mentioned above, the reader uninterested in this 
property may ignore the mention. In the critical step where the 

tra I sfor ? a I aon is inescapable, the galilean invariance of 
the pertinent variables is essential to the analysis. This step 

appears in the section entitled Geometric Acoustics and Blokhint- 
sev Invariance. 
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wavfi fronts and phase.- According to the basic concepts^of 
■■ t". 1 ■ “ — j- hp q i ffn al is 'propagated on wave fronts. Wave iront 

variable, as it will serve equally well to parametrize unew 

Sonfs’ what° it dimen si on s ? may Se ^fuhmp^Jnt. A phase may be 
chSse^to be dimensionless _ or to have the dimensions of time or 
distance, as may be convenient. 

As defined, the phase parametrizes the wave fronts^over the 
entire history of the wave propagation an , this global 

phone xixe my .,,0 time measured from the passage of a 

wave system goes by. Thus, time measui «u turns out to be 

reference wave front is a phase variable, one that turns oul 

efobal as well as local in a steady atmosphere; this particular 
glooai as _weij- „ hfll1 1]se in our general treatment of 

variable is the pittance measured at a given instant normal 

geometric acoustics . J^t erence wave front is a suitable local 
to the wave fronts fr re ference Mach cone in a coordinate 

phase. Distance , tn the aircraft is another local phase, 

ZlTs thfon^wS shSfuse in treating the flow near the aircraft. 

mo ill w pirate the distinction between phase and a local phase, 

r s n ?L or frr^sre nl Tw^ 

H S ^a-fs^rivriho^r^e 

used as a local phase. 

The reason the concept of phase is important is that we m ^st 
correctly id ®ntify the 2^i e ° V ^rSctly defined globally, 

serve^precisely ' thts^S^of^In our presentation of the theore- 
tical analysis, we use two local phases (L/L a and ) 
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one We raa y note that in the general case In 

which the atmospheric properties change with time, it Is Impossible 
to use a physically defined entity (e.g., time measured by a fixed 
observer) as the. (global) phase. In this case the phase must be 

physicai a interpretation? ^ ° TO Wlth valid 


Dhaqp W a particular reference wave front as the front of zero 

? h : wave front is a surface in space which is tangent at 

the aircraft to the Mach cone with vertex at some specified 
reference. point on the aircraft axis. Such a surface is termed a 
hP of C ?h° ld ' Po g convenience, we consider the reference point to 
nose ° f the aircraft. other points on the aircraft axis 
nhaq f vertices of Mach conoids or wave fronts of different 

5S rnJmf a C ^ 6P ' ts are discussed further in the section on 
i iacn Conoids and Gnound Intersections . 

Geom etric acoustics and rays .- The basic concept of geometric 
acousITcs.is that the signal is propagated along rays Ra|I are 

f?ont Ct and e ?b° f P ° int ? m ° vanS ln space • Each ray moves with a wave 
front, and the concept of the propagation of a signal on ravs is 

with that of its propagation on wave fronts. Since a 

7a Sn?n?°n?fV? i)e thafc ls ' a specification of the motion 

iLt-S tV Wlth d irae ^ at 1S a kinematic rather than a geometric 

P r r- appears desirable to emphasize this character, we 
a ra f a kinematic ray. The path of a ray is a geometric 

np-f-h ^ number ra y s traverse the same path, we term the 

flnrGpp geometric ray Since phase is constant on each wave front 
a lL' eac I ray moves with a. wave front, phase is also constant on 
f y ‘ . , , S en( ~ral solution the rays form a three -parameter family 

of point trajectories. The three parameters are analogous to 
L^grangian coordinates for particles moving in a fluid flow. One 
? he para meters is the phase, while the other two are selected to 
be an azimuth angle 0 and a time t a (to be defined later). 

,, In general, the rays corresponding to values of the phase 
other than zero. do not follow the same paths through space as do 

that r ^ S ir b?ch W ?i Ch the phasa ls zero * An important special case is 
that m which the ray geometry is steady, in which every ray path 

1S th ® path for a one -parameter family of kinematic rays. In this 
case the ray paths are what we have termed geometric rays which 
form a two-parameter family of curves in space. In applying the 
ana ogy to particles . moving in a fluid flow to this case, the flow 
is assumed steady, with the geometric rays then analogous to 
streamlines. The. property of steady ray geometry is not galilean 
invariant, and this fact indicates that the assumption of this 
property must be made with care. 


been 


Historically, for the most part only this special case has 
considered. Moving sources are rarely considered in geometric 
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oDtics- moving acoustic sources are generally treated in a _ coord 
inate frame in which they are fixed and, generally, only aircraf 
in steady flight have been considered as generators of sonic oom. 
Sul, historically, rays have been considered primarrly as geo- 
metric entities. 

We make the assumption of steady ray geometry in the sonic 
boom problem, with t a and 0 as the P^ameters for the geometric 
rays. This assumption is justified by the thinness of the enti 
wave system of interest, essentially by the fact that the.aircralt 
length is small compared with other macroscopic characteris 1 
noaies A ray emanating from the tail is simply so close to the 
corresponding one of zero phase that the dlffe ^/Hcknesfof 

the^wave^ system^and. te R* aHcr^scopic scale -asure the required 
condition is L a /R « 1 . If X ^ a measure. °f 
istic wave length of the acoustic signal, it is the condition 
X/R < < 1 which justifies geometric acoustics. The sonic boom 
problem is unique among acoustic problems in having X »L a with 
the consequence that the steady-ray-geometry assumption is valid 
whin gelme?r?c acoustics Is valid. Thus this assumption is sound 
even though the problem with a maneuvering aircraft is not a 
steady one. 

This assumption is basic to our analysis. It permits our 
calculating only the two-parameter family of rays corresponding to 
ztro Dhase considering the aircraft to be a single moving point m 
space Another basic assumption is that the cumulative nonline r 
effects do not affect the ray geometry. This is discussed later 
when we treat the nonlinear distortion. Hence, ray calculation 
follow linear theory. 

Besides the concepts of wave fronts, phase, and rays, another 
basic concept in geometric acoustics is that of ray tubes and ray- 
tube areas. Although ray _ tubes may be ^fined in the -general case, 
thev are much easier to visualize with the steady ^y g Y 

assumption. In the neighborhood of a given geometric ray we 
visualize a tube of geometric rays, i.e , a ray tube. The 
ponding; entity in the analogous steady fluid flow is a streamtupe. 
Hay-tube area is a measure of the differential area intercepted 
bv a surface cutting the ray tube and may be considered a vecto 
ISaltllv Lile a streamtubi, a ray tube is a differential quantity, 
SI I rSy-tubearea is actually defined in terms of derivatives 
with respect to the ray parameters. 

An element that greatly simplifies the calculation is the 
assumption that the atmosphere with its Hof°SnSil^s law 

£ SlStSic Swls then holds. This law permits the calculation 
of both the rays and corresponding ray-tube areas by quadratures. 
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hoc — - p St °5 the theoretlc a l analysis .- We turn now to a general 
description of our analysis of sonic boom, the details of which are 
presented m the subsequent sections of this chapter. The analysis 
be conceptually divided into three main parts which will appear 
in sections of the chapter preceded by a short section on The 

Dar? S nf e ih' at ^ he . end we add a Note on Viscous Effects. The first 
p.art of the analysis comprises the sections entitled Aircraft 

sections S \S^ aal r WaVe JJ ormals ' Mach Conoids and Ground Inter- 

cerns the Sinl W a ? /-? ay Tracln S; aa d Ray-Tube Area. It con- 
cerns the calculation of the rays and ray-tube areas for zero phase 

entlt^d e p? nCe M PhaS ^* The second part comprises the sections 
entitied Fiow Nesir the Aircraft, and Geometric Acoustics and Blok- 
hintsev Invariance. It concerns the calculation by linear theory 
o acoustic signals along each geometric ray. The third part 

222^“ Location 10 "?^"" 1 " 16 " \ 1Snal Portion anflge "variable , 
and Shock Location. It concerns the calculation, with shocks 

properly accounted for, of the nonlinear distortion of the signal 

Otvector quantities are introduced and used in thl ^ 

analysis The components of these vectors in the various coordi- 

Components mS USed Er6 glven ln the sect ion entitled SYMBOLS, Vector 

The maneuver of the aircraft (strictly speaking of the w fpr> 
ence point) is required in detail/ Variable s ^fjnt reduced In the 
snapper A1 ^raft Maneuvers which describe the trajectory i^ 

the orientation of the flight axis, the velocity of the air- 

all L^uS atmo ®P here ' and ^e local sound speed, 

™nnM fUnCtl ° n n S ° f tlme ta * Tlme derivatives of certain of the 

calculation"^ f° r later USe ln the ra N- tube area 

ulation. At each instant t a we visualize a Mach cone attached 

to the nose of the aircraft. Thl normals to the Mach cone form f 

e parameter family of directions forming a wave-normal cone with 

Si T a ^c e L beinS an azl r th angle * • The two quantities u 
and 0 are the ray parameters discussed earlier. 

tde sec tion Mach Conoids and Ground Intersections, the wave 
ronts and rays from an aircraft in maneuvering flight are dis- 
cussed generally, with particular attention to the interactions of 
the rays and wave fronts with the ground. intersections ol 

• • +. . generators of the wave-normal cone at the aircraft are fhp 
initial wave normals for the calculation of the rOyO ?he oSL- 
tation of these normals is known as a function of the ray para- 

i F ° r eaCh WaVe norraal we calculate two quantities which are 
invariant on rays according to the appropriate Snell's law These 
invariants are then used to calculate the ray trajectori^.' 

cuttinff e n?-?^ Ube A ar ! a ds deflned as tha t given by horizontal 
cutting planes . An analytic expression for this area is obtained 
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in terms of the maneuver variables, certain of their time deriva- 
tives and three quadratures along the ray. The ray- tube area is 
£hufib“ined as a function of altitude along each ray and may be 
calculated concurrently with the ray trajectory. 

For the second part of the analysis (Flow Near the Aircraft^ 
Geometric Acoustics and Blokhintsev Invariance), consider fir_^ 

the flow close to the aircraft. In particular, we need the asymp 
totic form of the local solution, valid at a distance from the 
T?^ht axis large compared with the effective lateral dimensions 
of the aircraft but small compared with characteristic scales for 

the atmosphere. This asymptotic form of the d ° c ' al A S ^ ^ff iciently 
interoretable as a geometric acoustics solution. At a sullicien y 
^ar-e distance r in a particular direction away from the flight 
axis of the aircraft, the solution appears the same as that from 
Sne distribution of sources and sinks, the same as that from an 
equivalent body of revolution The pressur^ perturbation in th 
asymptotic solution is proportional to r V times a function 
of a suitable defined phase and of an asimuth angle 0 r (simply 
related to 0). This F- function depends also upon the Mach num 
and lift coefficient of the aircraft, which are functions of th 
time t The F-function is then a function of phase and of 

ray partmhers t a and 0 and is invariant al °^ t ®f^ c , k ““ a h C 
ray It is obtainable either by a computation (outlined here in 
the section Flow Near the Aircraft) or from experiment. It 
assumed to be a known function in the computer program. 

Tn the general stratified atmosphere with winds, the appro- 
nriate genlrtl definition of phase is as the time £ measured by 
an observer fixed in a ground-based coordinate system, and define 
to be zero the instant the zero-phase wave front passes. Invar 
iance results of Blokhintsev permit the acoustic signal of ea 
ray to be described in terms of a function V e (£) which 
constant on the ray. 

The relation between the local and general phase variables is 
then found so that F may be expressed as a function of 4 The 

relation between the function F f 0 ^tJ^rav (ti£ce F is assumed 

known) 1Ve The h relation°betwien Y e and pressure perturbation Ap 
is known, so that Ap(£) is determined at any point on each 
geometric ray . 

Tn the third part of the analysis (Signal Distortion and Age 

Qh-ittc; in the signal, whereby a given point m the signature 
Sy^PPeaflaSiS than predicted by the linear theory. 
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d i stort ion arises because compression waves travel 
slightly faster, and expansion waves slightly slower, than do 
in ^ mitesimal ^ disturbances . In terms of the phase variable £ 
this phase shift equals V e times an age variable t which can' 
be computed along each ray by a quadrature. The distorted signal 

Wonafto 3 x ° rlglnal ° ne V e(^) sheared by an amount propS- 

severa^ e v ein^ r ^ d ^ lgnal may be rault i va lred and may thus give 
several values of the pressure perturbation for a single value of 

£ . Ihysically, this anomaly indicates the presence of shock 

A leoaSte skocks are Properly taken into account. 

an ®; 1 y ais sh °ws where shock waves must lie and shows 
™ bach P f rts of the signature have been "eaten up" by the shocks 

t±Ll°e aPPear ‘ The reSUlt ° f the anal ^is is the complete, 
Shock^ d -S r eu SUre sl S nature at any desired point, with 

tu if they are present. The nonlinear effect does 

the P ressure perturbation insofar as the 
longer °app ear ? ri ^ ina ^ Slgnal th&t are eaten Up by shocks no 

the rav e tnhp°JLf a i 1S near & caustlc ; a surface in space at which 
the ray-tube area. becomes zero. It also fails near the boundary 

a shadow zone into which no rays penetrate and may fail near a 

f°r. WhlCh the function is singular in some way. 

The lmes-r solutions m these regions are solutions to diffraction 
problems for which geometric acoustics is invalid These problems 
are outside the scope of the analysis of this report! P 

The Atmosphere 

. The coordinate system used is cartesian with x and v 
onzontal distances east and north, respectively, from a refer- 
^ 2 ? on J he ground . The ground is assumed level at an 

£rron£d de afc ove sea level, and z is altitude above the 

ground. The atmosphere is assumed to be a calorically perfect gas 
(constant specific heat ratio 7 ) with the thermodynamic proper- 

of 3D^d Pe » T ia# nslty P ’ P ress “re p = RTp , andhpelS 
of sound a = (7 e RT)V^ g lven as functions of altitude The 
pressure obeys the hydrostatic law " * 


dp 

dz 


■pg 


( 1 ) 


f the acceleration due to gravity. Winds are horizontal 
magnitude . u and direction dependent only upon z . The 
wind direction is specified by the wind heading angle r" measured 

Mnn f -°I? S° rbh - In accord wlth a ^±ent historical cSn?ln 

tion, the wind heading is taken as the direction from which the 
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wind comes ("the 
this convention, 
ponents ( _u x^ _u y 


north wind 
The veloci 
) with 


doth blow"), and we have acceded to 
ty of the wind has east and north c 


om- 


u = u sin r] 
x 


u = u cos T) 

«y 

. r m 1 cation in the program .- Inputs into the ^gram include 
thP t emperatu re — T , either the density p or the pressure p , 
the wind speed u , and the wind heading gl ^®g 0 lnpu t is 

functions of the altitude Z+ Q provides for specification of 
Se §6^.1 Stafdard G Atmosphe?e (re?. 2) with any wind distri- 

bution . 

There is no provision in the program to ensure that the hydro- 
statifrelation (?) is =^isfied Hydrostatre consistency of the 
input data is the responsibility of the operator. 


Aircraft Maneuvers 

Tn our study of sonic boom, we need the trajectory of the air- 
craft in order to know where the rays i^not^pecif led , are 

integrated for the aircraft ^^“^Ai^atives of heading angle, 

f if anfaf craft ' velocity or J e 5St???.g e ?he ater 

calculation"©^ the se^ derivat ive^ from KM load fetors are 
also included here. 

The aircraft moves through space supersonically on some known 

x a ( t a ) > ya(t a Ji ana nose. The subscript a 

identifies^var tables 6 defining the aircraft position in a ground- 
fixed coordinate system. 

The first stage analysis is^n ^stigation^of ^h^^ 

equations governing the traj^ J^ lnate syst em in which the local 

hJB velocity L ™est U e , in a coordinate system moving with the 

Ihis veiohty has magnitude V a head ng angle 

°Thf direction Iff’ If a^nSanths tlrmed the 

flight axis. 

With respect to the ground-fixed coordinate system, the diff- 
e rent ialfquat ions for the flight trajectory are then 
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WT = V cos 7 sin f - u (z ) 
a a 

dy a 

ItT = V cos 7 cos 0 - u y (z a ) ( 2 ) 


dz 


a 



= V sin 7 


If 

the 

y a 


V ^7 , and f are known 
third equation to obtain 
by integrating the other 


functions of t a , we can 

z a(t a )^ and then obtain 
two . 


integrate 
x a and 


The 
divided 
path is 
defined 


acceleration of the aircraft is equal to the net force 
y he aircraft mass W/g. The component along the flight 
b <, n T - sm y)e, where n T is a net thrust loaf fiftor 


n 


ij, — ( T D)/W = ( C T - C 


D 


^co S ref/ W 


(3) 


D representing the thrust and drag on the aircraft 

I ctricl C Ta thf^ f.TrPnfff QnH hq TO CT /nr» /r ^ -f 4 -P 4 -! — Jf _ _ _ _z_ 

Qco 


with T and 

respectively, Crp and the thrust and drag coefficients 

the dynamic pressure (l/g)p V 2 , and S ref thf ' 

»^r¥e- are i T he entity g cos 7 r fl a component o? ?he 
fligldlflf. graVlty a ° tlnS laterall y' normal to the 

side 1 flrees 01 ’and to he^a®? !° 5 e laterall y symmetric, without 
axis 1 in?? d t ^ banked at an angle 0 a about the flight 

t-iiS' ThS ll £ t ° n the aircra -ft then provides a normal accelera- 
tion component n L g , where n L is the lift load facto? de??ned 


n. 


= L/W 


C L q oo S ref // ' A? 


w 


L mt^ nd - ° L . re P r ' ese nt the lift and lift coefficient respect- 
y. This is directed so that n^g cos 0 a opposes the gravity 

showTth °? S h- wlaile n L§ sin 0a is horizontal. Figure 1 
shows the acceleration components while figure 2 shows how^ rf> is 

inatfframe ^ows the lateral acceleration components. The coord- 

i^lativp fn thi iL' ] ls l otatad */2 - 0 counterclockwise 
relative to the reference frame (x,y,z) and is used laten to 

develop wave propagation directions. 
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To obtain 
components , we 
combine terms. 


expressions for the lateral 
differentiate equations (2) 
This yields 


and axial acceleration 
with respect to t a an 


,2 

a x 


n T g sin 0, 

J_j c 


COS f p 

dtjf 


a 


- sin "0 


,2 

d y 


dt 


a 

2 

a 


d^z a 

(n T cos 0 - cos y)g = cos 7 o ” 

' L a dt a 


,2 

d x 


a 


2 



sin 7 


sin ^ 


+ cos 0 



Figure 2. _View looking forward along 
acceleration components, bank angle. 


flight axis showing 
and azimuth angle. 


( n T - sin 7 )g = sin 7 


d 2 z 


a 


dt 


2 

S 


+ cos 7I sin f 


d 2 x 

Si u 

2~ + COS Ip -rf- 


<3 2 y„ 


dt 


d t 


a 


Carrying out the differentiation with respect to t 
obtain a 


, we then 



V cos 7 


dz / du du \ 

§|- = n L s sln 7 - Tfr( sin * di 2 - oos * d5~) 

a. a 


V 


d7 _ 
dt ~ 


(n T cos 0 - cos •y] 

1 i cx 


( 5 : 


dz 


cl 


d t 

cl 


du 

sin y( cos ip + sin f 


du 


x 


dz 


dV 

dt 


dz 


(n T - sin 7 )s + cos T^cos ip + sin ip 


a 


du 


du 


The factor dz a /dt a 7 rorn ( 2 )j simply V sin 7 . Equations 

(5) relate the two load factors n L and n T with the time deri- 
vatives of ip , 7 3 and V • 


Although they are not directly involved in the aircraft dynam 
ic s , the speed of sound a and the aircraft Mach number M = V/a 
are convenient to use in the acoustic analyses. The speed of sound 
is obtained from 


a 


( 7 e R T ) 1 / 2 


( 6 ) 


The gradient of the speed of sound satisfies the relation 


da = ]e^ dT 
dz 2a dz 


The time derivative of V = Ma may then be expressed as 


dV 

dt 

cl 


dM 
d t 

a 


+ 


dz 


M 


a 


dt. 


a 


da 

dz 


= ai 


( dM 
V dt 


+ M sin 


•y da 
y dz 


(7) 
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Application in the program .- Inputs into the program include 
and M as functions of the time t . Equations (2) are 
integrated to obtain the aircraft trajectory. The time derivatives 
of f , y > and M . (actually of the Mach angle p. = sin.- 1 (l/M) 
are required later in the program for the ray-tube area calculation 
and,, for this purpose, two Maneuver Options are provided. In the 
first option., the load factors n L and n T are additional inputs, 
and the time derivatives are calculated from equations ( 5 ) and ( 7 ). 
In the second option, the time derivatives are calculated directly 
within the program by differentiating the input data. 

The purpose of including Maneuver Option 1 is to provide a 
more accurate calculation of the time derivatives in case the load 
factors are accurately known, as perhaps from accelerometer data 
from . a flight test. In this option, the input data are redundant, 
and it is the responsibility of the operator to ensure that they 
are reasonably consistent. 


Initial Wave Normals 

The purpose of this section is to express the initial orienta- 
tion of the wave fronts as they leave the aircraft. This initial 
orientation gives the basic parameters needed for the ray calcu- 
lation. 

Here we have an example of the principle discussed earlier of 
using galilean invariance to identify the variables which are 
preferable for use in the analysis. The wave normals are the ray 
directions in one particular coordinate system, that fixed in the 
undisturbed atmosphere at the aircraft altitude. Ray directions 
are not galilean invariant, while wave front shapes and wave 
normals are. The use of wave normals rather than ray directions to 
define the basic variables keeps the analysis in its simplest form. 

^ ith each instant of time t a during the aircraft flight, we 
associate a Mach cone with vertex located on a given reference 
point on the aircraft. This Mach cone is tangent to the Mach 
conoid (wave front) moving with the aircraft. The normals to the 
Mach cone at the vertex form a wave normal cone. We consider an 
instantaneous coordinate system (x^y^z) for purposes of descri- 
bing direction only (fig. 1 ), rotated an angle tt/2 - 1 ]/ counter- 
clockwise relative to the basic coordinate system, with origin at 
the reference point on the aircraft. The two cones are illustrated 
in this coordinate system in figure 3. The half -angle of the wave 
normal cone is the complement of the Mach angle q. . 

An arbitrary azimuth angle 0 is chosen, measured according 
to a right-hand rule about the flight axis from the downward 
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V--P 





direction (see figs. 2 and 3) • The particular wave normal corres- 
ponding to this azimuth angle is considered and will determine the 
particular ray determined by the parameters tg^ 4 • The direction 
cosines (£,m ,n) of the wave normal relative to the axes (x' y' -z) 
are to be determined next. J 

The pertinent vector directions may be represented by points 
on a unit sphere. The desired direction cosines can be derived 
either by spherical trigonometry or by calculating cartesian coord-* 
inates of points on the unit sphere. Choosing the first approach 
we consider the spherical triangles of figure 4 and obtain " ' 



Figure 4. Spherical triangles for calculating wave 
normal direction cosines. 
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£ = sin |jl cos 7 + cos \± sin 7 cos 0 
m = cos |jl sin 0 

n = - sin |jl sin 7 + cos \i cos 7 cos 0 


( 8 ) 


using the law of cosines to obtain f and n and the law of sines 
for a right triangle to obtain m • Essentially the same results 
are obtained in Appendix IV of reference 3 . The downward or (-z) 
axis is used here because we want to consider prxmarxly descendxng 
acoustic signals and rays. 


The angle v is defined as the heading angle of the wave 
normal (fig. 3(b)), and the angle 0 O is defined as the angle of. 
the wave normal below the horizontal (fig. O • Using these defxnx- 
tions, m /i = tan (f - v) and n - sin 0 O • With equatxons (8), 

these yield 


v 

and 


sin 0 

o 


The maneuver history of the aircraft provides p and 7 , so that 
these equations give v and 0 O as functions of the two ray para- 
meters t a and £ . We note that 

£ = cos 0 Q cos (f - v) 

m = cos 0 Q sin (^ - v) 


Another coordinate system (x-^yq^z) is shown in. figure 3 and 
is aligned with a particular wave normal. This coordinate system 
is not used in this section but is used below in the treatment of 
ray tracing and ray- tube areas. 

Application in the program .- Using the known values of f s 
7 , and M , the quantities v and sin 0 q are calculated from 
equations (9) and (10) as functions of the ray parameters t a 
and 0 . 


-if cos p sin 0 

f - tan "l. sTn p cos 7 + cos p sin 7 cos 0 


- sin p sin 7 + cos p cos 7 cos 


= n 


( 9 ) 

( 10 ) 
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Mach Conoids and Ground Intersections 

At this point we are ready to calculate the rays. The purpose 
of this section is to describe the rays and the Mach conoids (wave 
fronts) in general terms before going into the detailed calculation 
of the rays. The intent also is to show the functions needed to 
describe the rays and wave fronts globally and needed to determine 
when and where the sonic boom signals hit the ground. The primary 
purpose of the section is, thus, largely conceptual, and the reader 
primarily concerned with the algorithm may skip the section. 

As the aircraft moves through space, a wave system associated 
with the aircraft moves with it and propagates away from the flight 
path. The wave system, which consists of a one-parameter family of 
wave fronts, is characterized by a single wave front chosen here 
to be the one of zero phase. This wave front is attached to the 
aircraft at the reference point on the aircraft and is tangent to 
the Mach cone associated with this point. The front is the same as 
the Mach cone only in the special case of straight flight at 
constant speed in an atmosphere of uniform temperature. This 
reference wave front is termed the reference Mach conoid and is 
shown schematically in figure 5. 

The reference wave fronts or Mach conoids are not calculated 
directly. What we calculate are the kinematic rays corresponding 
to that. wave front, which here are the rays of zero phase. With 
the basic assumption discussed in the section General Description 
the ray paths for these rays are also those for other values of the 
phase and have been termed geometric rays. 

The rays are specified by three functions giving x , y and 
z as functions of t a , 0 , and t . Here t is the time on 
the ray, while. t a and 0 are the ray parameters defined in the 
preceding section. In a stratified atmosphere we replace t by 
z as. independent variable. The ray is then specified by the two 
functions x(t a ,0,z) and y(t a ,0,z) describing the ray path (or 
geometric ray), together with the function t(t a ,0,z) giving time 
along the ray. The calculation of these functions is described in 
the following. section. The ray paths for a given time t a of 
emission are illustrated schematically in figure shown here as 
straight lines. In the general case, they are curved lines. 

The ground is at z = 0 . The solid ground intersection curve 
in figure 5 is. given by the functions x(t a ,0,O) and y(t a ,0,O) 
for a. given emission time t a . The time of arrival of the signal 
on this ground intersection curve is not constant and is given bv 
the function t(t a ,0,O) . 

The wave fronts or Mach conoids are surfaces of constant t . 

An inversion of the function t(t a ,0,z) gives a function t a (0,z,t). 
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Figure 5. Maneuvering aircraft, showing coordinate system, 

rays, and Mach conoid. 


When this is substituted into the ray path functions, the result 
is two functions, x(0, z,t) and. y(0,z,t) , giving the Mach 
conoids parametrized by the variables 0 and z . 

The intersection of a Mach conoid with the ground is the 
curve given by the functions x(0,O,t) and y(0,O*t), obtaina le 
by the process just described but with z _kept equal to zero 
Such a ground intersection of coincident signals is illustrated y 
the dashed curve in figure 5- 

The phase variable which is introduced later in the analysis 
is a time £ measured from passage of the zero-phase wave front. 
Rays and wave fronts of phase different from zero obey the same 
equations as do those of zero phase, except that the arrival time 
t is replaced by the time t + £ . 


Application in the program .- After the rays are calculated, 
following the procedure described in the next section, .the Mach 
conoid ground intersections are obtained by interpolation o e 
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variable t . These ground intersection curves, together with 
corresponding maximum values of the pressure in the final signature 
are part of the output . 


Snell's Law and Ray Tracing 

The purpose of this section is to present the appropriate 
Snell's law for geometric acoustic propagation in a stratified 
atmosphere and to derive with this law the equations whose quadra- 
ture gives the rays. 

^ -^- n acoustic theory a signal on acoustic disturbance is propa- 

gated on a moving wave front. (Prom the mathematical point of 
view, a wave front is a characteristic hypersurface in four-dimen- 
sional space-time for the full equations of motion.) It moves in 
such a way that its normal velocity relative to the medium is the 
speed of sound. Its actual normal velocity in space is 

c n - a + n . u ( 11 ) 

where .n is a unit vector normal to the surface pointing in the 
direction of propagation and u is the vector velocity of the 
undisturbed medium (wind vector). 

A signal initiated at one instant from a point is found a 
short time 6t later_wit.hin a spherelet of radius a6t whose 
center is displaced u6t from the original point. If every point 
on a wave front emits a signal at a given instant, Huygen 1 s prin- 
ciple identifies one of the two envelopes of the spherelets 6t 
later as the wave front at that time (the other envelope corres- 
ponds to -n and is usually without meaning). This principle 
gives a motion to the wave front in accord with equation (ll). 

In geometric acoustics, the concept of a ray is fundamental. 

A ray is a point trajectory and may be defined 

(a) as a characteristic in an asymptotic development of the 
equations of motion for small wave length; 

(b) as a bicharacteristic for the full equations of motion, 
corresponding to the wave front as a characteristic 
hypersurface; 

(c) to . move from the point of emission of a spherelet to the 
point of tangency of the spherelet with the envelope 
wave front at a time 6t later. 

Any of these definitions leads to the result that the ray is a 
trajectory of a point that moves with the velocity 
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dr 

dt 


= c = an + u 


( 12 ) 


where r - xi + yj + zk is the position vector of the point and 
I T and k are cartesian unit vectors. This velocity is 
termed* the ray velocity or group velocity In geometric acoustics, 
the signal is propagated along rays with this velocity. Note that 

c n = H • 5 • 

Besides the ray or group velocity, an Important entity is the 
slowness vector or inverse phase velocity n / c n • # Thi y 

more familiar in geometric optics than in geometric acoustics, 
primarily because, the subject of geometric optics has been so 
thoroughly studied and applied to practical problems. 

In order to calculate rays, it is necessary to know how the 
wave normal vector £ changes along rays A general refraction 
law may be derived (see ref. 4, for example) which states that 


dn 

dt 


Va + (Vu) • n - n[n • Va + n • (Vu) • n] 


(13) 


alone a ray. The combination of equations (12) and (13) i ® a . 
system of differential equations which must be solved to obtain 

the ray . 

In the case treated herein of a steady horizontally strati- 
fied atmosphere, the calculation of a ray is much simpler A part- 
icular refraction law or Snell's law is available which g 
pirn licit lv alone a ray (see ref. 5; for example). This Snell s 
law", staged in Its mit general form is that the horizontal vector 
component of the inverse phase velocity vector n/c n is constan 
along each ray. We decompose n into horizontal and vertica 
components according to 


n = cos 0 n ' - sin 9 k 


is a horizontal unit vector. The angle 9 is the 
n below the horizontal. Our Snell's law then^state 


s 


where n ! 

angle of n below the horizontal uux- oncxx-o ~~ 

that the horizontal vector cos 6n/c n is constant along eac. r y 

r -7 \ T.T^v ‘ ~ ”* ° 


that the horizontal vector cos u n/c^ is oe 
(see figs. 6 and 7). We define the velocity 


as 


'n 


t -\ ) i \ 


in terms 
initial 


of which the invariant horizontal vector is 
value of 9 when the ray is emitted from the 


- 1 - . 

c© n 

aircraft 


The 

is 
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Figure 7. Oblique view of velocity plot. 


the ano-le 0 O defined by equation (10). The quantity 


, and the 

direction of n' arrtiSsiOTah'anralong each nay The direction 
direction dpnoted v in accord with the notation 

or heading angle of n is denorea v ±u » 

of the previous section. The two invariants c Q and 


1UU O VV o/ 4. X V v^v J- * w . 

v(t a ,0) Of the two ray parameters. 


functions c o (t a >0) and v(t a 

We introduce a coordinate system ( x l>yi<7^) ?ii Sn ??hU that 
the wave normal n lies in the (x^-z) plane ( 

relativ^to tL°SL?f coSrdSa^systlm". ' The main -e of this coor- 
dinate system will be in the following section, in the calculation 
of ray-tube area. 

The wind vector u has components (-u x ,~0 ln ( x ^y) 

frame and components (-u n ,-u t ) in the (x^) frame with 


= u cos (v - T\) = u x sin v + u y cos v 
= u sin (v - r]) = -u x cos v + u^ sin v 


(15) 


The minus signs before the components come 
convention mentioned earlier. Note that 


_from the ancient 
u • n = -u cos 6 


wind 
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For each initial wave normal (and corresoondi up- Mine <~>-p ^ \ 
m the wave normal cone at the aircraft th™ S of 0 ) 

ray. The value of v for this Sv ?, hJ ?' ILl °"®,°orrespondlng 
(9) ¥e obtain n -r^ +-v, niS r ~ y 1S that obtained from equation 

{?). we obtain c Q for the ray from the expression 


'n 


o 


c - 

0 COS Gr 


a + u 
o o 


n 


o 


cos 0 


o 


a 

o 

cos 6 q Un o 


( 16 ) 


where 


o 


re cos e Q is obtained through equation (10) and the suhscrint 

s sS~- -« 


COS 6 



(17) 


ffSncMon^T f and M r °V°° - and = and thereby as 

bp f wppn , n a t ’ and z • figures 6 and 7 show the relation 
m b * > and a and the wind components -u n and ~u 

The wave front appears edge-on in the (xt z) ni nt +- -u * 

b;fhas m no n J?feot" U h th^sSen's^aw^ fr ° nt d06S affec * 5 6 


we 


In 

need 




(see figs. 6 and 7) : 



dy 3 

dtT 


-U 


t 


d(-z) 

dt 


= a sin 9 


These equations are components of the vector equation 1 12 ) The 
equations"^ VaPlable 1S Chan£ed frora * to * ^Aeke ^ 


dx^ 

d(-z) 


a cos 6 


u 


n 


a sin 6 


dy- 


-u 


t 


d( -z) a sin 6 
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A transformation 
equations 


from (x 1<yi ) to (x,y) gives the desired ray 

dx a cos 9 sin v - u x 

d( -z ) = a sin 9 


dy 

d(-z) 


a cos 9 cos v - u y 
a sin 9 


(18) 


dt = 1 

d(-z) ~~ a sin 0 


With a downward propagating ray (sin 6 > ^^^thrini- 

equations (l8) is carried out xn the -z directxon irom 

. _ „ . f ~ ^ v _ ~\T — ^ 


t = t a ) to the ground 
are obtained as 


tial point (z = z a , x = x a , y - ^a > 

(z = 0). The functions sin 9 and cos 9 
functions of z from equation (17)» 

Foliations (l8) are the basic equations of this section _ as they 
functions IfTTZT) XlTZltT TiSculsed “the* 

modification to take neighboring rays into account. The ray 
eauations in this form are more convenient than equations (18) fo 
computing rSy-tibe areas and will be used for this purpose. 

Application in the program .- In the program, for the selected 

values of the ray parameters (ta>8) , c o is calcula e Eaua tions 

equation (l6) with u nQ obtained from equation (15) ^ E ?^ lons 
('18') are then integrated for the rays with 9 obtained iiom 
innit inn ( 17 ) Only downward propagating rays are calculated, 

£e calSulahin £ stopped when the ray is approximately horizontal 

Historical not| - The faction law of Rajlelg £ fs ^ 1 ' 8 

re?! 6) [nk.Tr fill. Rayleigh £d not distinguish 

between wave normals and rays, however. . J^ton, W as n?t i£ thA wive 
4-Hnt the correct ray propagation velocity was nor in one wo. v 

normal^direction? r He gavl the correct planar ray tracing equations 

T.ri tt PYamnles Fuiiwhara, m 1912 (ref. o), gave tne corr 

Snell’s law iA three dimensions and the J or ^f s P°^ntif ie^thr^y* in 
eauations (18) with examples. Emden ref. 9) identified the ray 

terms of energy transport without defining the energy. 
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Ray -Tube Area 


rav-tube area°J?on f sectlon ls to obtain an expression for 

ray tube. area along the rays. The ray-tube area is needed subse- 
quently in the analysis in order to express the acoustic signal 

di“ort?on! ely ’ and alS ° thereby t0 ca?culate thS nonwLS 8 al 

tube made y un U of if & d jf f ® rentlal concept, to be visualized as a 
ruoe made up of rays which are differentially close to the nar=H- 

suih r a r tSbe e js S th nVeS ^ 1 Sate H' The actual cross -sectional area of 
, ? h , a tude 13 thus also a differential quantity The quantity 

iS a finlte ™ea S ure y h sSch^ «ffL- 

a d n0t actuall y a Physically identifiable area. 

di ffpp lyinS - an L ra r tube area by a constant factor will make no 

this reSUltS ° f the anal y sls - A consequence of 

his fact is that the dimensions assigned to the ray-tube area are 

compieteiy unimportant and may be changed through such a factor to 
mu^ipMca??In n ^ n ?h' The invariance of the final result to a ° 
was used as a check SflL anSlysIt.^ “ arbltrap y instant factor 


In our case, we use the 


duced earlier, correJponikg ^h^’p^ticSJafSngle^v^h 1 ^?^ 

parameters 110 ! th n r&yS P aramet; nized in terms of the ray 

parameters t a and 0 . Our use of Snell's law for a stratified 

mosphere directs the use of horizontal cutting planes with the 
vector ray-tube area directed in the direction of the -"Sis We 

area V on U a 1 DlanP dlfferent j al nay-tube area as the quadrilateral 
meter° ft T) f = coa f an f determined by the rays with para- 

« “uuB^i it\ + 60) 

de«L t ^ J ^ y 0b -iuSe°a f rei X h y h s (W) f^Se 


A(t ,0,z) = 


o 


dx-j^ 


^a 

^a 

dx^ 



<^0 


(19) 


in terms of the functions x i(t a ,0,z) and y-|(t~,0,z) These we 

fSnctions 1V xft f f? a r °tation of - V p from the 

• -u X (^aj and y(t a; ,0 J ,z) that were obtained from 

fSrm Sra The S fl?tS i °c-l (1 - ) -- T he analysis is much simpler in this 
lorm. lhe factor c 0 i is included in the definition of A so as 
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nimrafi flight path 


+ 6t 


nt ial 
e area 

z = constant 


eighboring ray 
t + 6t , 0 + 60) 


Reference ray Neighboring ray 
(t ,0) (t a ,0 + 60) 


Figure 


Sketch showing ray-tube area. 


to make the subsequent f ^mula^simpler; ^the factor^ & here 
dimens ions^of ISnlth in this de 

Son^hfSf Sd eventuall^arrive^a^the emission ( S 6) below. 

A few authors have defined entities analogous £ ay ;£ u ^ e ®£ • t 

sss4s 
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fn™ W o a ^ m °^ e general terms, may also be found in reference 4 
inn^ J^ 1 ? 13113 ^ ay be consi( 3ered a special case of a general 
Tmw^° aCh t0 ° alculatln g Jacobians of analogous type for general 
being^prepared?) 1 ^ * ^ ° n th±S S ^ ect by the^rft^^L 

as in^a^LS 4° 

~ - &^ss.ts h ^ i £'~ S >* 

til Tn^^iedVX) where ° ntal Ve ° t ° r 3f ' /3t a (sh °” S “ 

r = xi + yj = x^n' + y,k x n' 


We rewrite equation (19) in the 


form 


c A 
o 


(^ii) 

Wo 

Cia 

Wo, 

+ 

St a 

(ti) 

ydt a y 0 

dy i . 

3t a W a / 0 

1 — 1 
X 
ro 

1 — 1 
!>a 

SO 

3x l 


^1 

1 M 

c>0 1 


b(p 


bep 


KS,r:,;:rrs;: s s s“sss.“s»*“,s;,;- 

“° e ?aJJ e r UPPe fH tern ? in the aefe^inant^e'inte'gSaJrover 

.•„g, nom^ the point of emission. We make a transformation of 


. - WJ. CUUDD, 

the independent variables (t a , 0 ) 


(c Q ,v) and write 


c A 
o 


0 

dv 

d0 


o 


/bx- 

^\bt 

bx. 

dc„ 


a'O 


dc dx- 


d0 dv 


by 1 

be 

bv 


o 


be N 

0 dv 


Sy 1 



dc o Sc o 

dC -y 

0 dv 


Sy x 

d0 d0 


dv dv 

- 


(20) 
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Here we have used the theorem that the determinant of the product 
of two matrices is the product of their determinants. The deriva- 
tives of Xi and y^ with respect to Cq and v in equation 
( 20 ) are again integrals over -z taken from the point of emission. 
Our next step is to evaluate these quantities. 


The initial values of the t a derivatives may be calculated 
with the aid of the equations for the direction cosines of the 
initial wave normals. We note that 

S, = cos 9 q cos (if - v) 

m = cos 9 sin (f - v ) 
o 

The result of the calculation is 


dx-A 

StJ 


= V[ cos 7 cos(^ - v) + sin 7 cos 9 ] 


a y 0 


u 


n 


o 


1 + 


sin 7 


sm |jl sin 


o' 




sin 7 

c o v + sin p sin 0 




~dt 


a/ 0 


/ , \ / sin 7 

V cos 7 sm (if - v) - u t ^ ^1 + s±n ^ sin-3- 


o 


V cos 7 cos \i sm 
cos 


o 


- u+. 1 + — T 


sin 7 


sm |i sm 


o- 


( 21 ) 


( 22 ) 


To calculate the derivatives of (xp,yi) relative to c G and 
v , we first recognize that the coordinate system is defined to 
correspond to one reference ray with v = and that we must 

consider neighboring rays. The ray tracing equations (see equa- 
tions preceding ( 18 )) are written in terms of X]_ and yq and are 


dx-, Q u 

1 _ cos 9 _ n 

d ( -z ) _ sin 6 a sin 9 


dy x 

d’(-z) 


cos 9 
sin 5 



u, 


a sin 
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In the second equation, (v - v r ) represents sin(v - v r ) with 
v - v r small. A factor cos(v - v r ) in the first equation has 
been set equal to 1. In the same terms, equation (17) relating 
c 0 and 0 may he written ” 


c 

o 


a 

cos 0 


(23) 


These three equations are differentiated at constant z with 
respect to the three variables c Q , v , and 0 , and then v is 
set equal to v r . The differential of 0 is eliminated to yield 
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We define the following integrals 


V-z) 


f cos^ 0 d (-z ) 

^ -z sin^ 0 

o 


I 2 (-z) 


u^ cos^ 0 d(-z) 

O O 

■ z a sin J 0 

o 


( 24 ) 


I 3 (-z) 


2 9 

0 zu, cos J 0 n \ 

/ / t , £0^0 4/ 

1 V 2 . 3 fl + sin 0/ ' ' 

-z s a sm 0 

o 


as indefinite integrals equal to zero at z - z Q . This notation, 

avoiding the dummy variables of definite integrals, is a conven- 
ient one. In terms of these integrals, we can evaluate the 
quantities 


JL 

c 


dx n 


dc 

o o 


= - I- 


_1_ 

c 


o 


dx 

^v" 


1 


dy- 


o 


dy n 

± — _ T 

x 3 


( 25 ) 
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We now substitute the expressions obtained in equations ( 21 ), 
(22), and ( 25 ) into formula ( 20 ) for the ray -tube area. The result 
is 


A 


1 + 


sm 


X 


sin p sin 9 


O' 


(I Q - I, ) - Y sln co s M- cos 7 I 

^ ^ o 1 cos 0~ 1. 


l 

bcj) 


1 + 


sin 7 


sin \i sin 0 


o 


(I 3 ■ u t 0 I 2> 


V sin 0 cos j j, cos 7 

cos 0 Q 2„ 


bv 

d0 


be 


o 


bv 


be 


o 


bv 


dt a b<f> bef) bt 


a* 


~ I2) 


( 26 ) 


This is the desired expression for A . This expression is given 
in terms of the derivatives of c Q and v with respect to the ray 
parameters, and these derivatives must now be calculated. 

In the calculation of these derivatives, the quantities i {/ , 

7 , and . p are functions of t a alone. The derivatives of v 

are obtained with some algebraic manipulation from equation (Q) and 
are 


6v _ diA sin 0 


dt 


dt 


a 


cos 2 0 


cos 7 ^ + cos p sin 0, 


o 


° dt. 


6v 

60 


( 27 ) 


cos p 
COS 2 0 


(cos p sin 7 + sin p cos 7 cos 0) 


To obtain the derivatives of c Q , we first differentiate equation 
(23) again at the aircraft (with z variable) to obtain 


a sin 9 


dc 


cos 2 0 


d0„ + u-j- dv + 


o 


da 


o 


cos dz 
o 


o 


du 


■n 


o 


dz 


dz 


Equation (10) for sin 0 Q is also differentiated, and the 0 O 
derivatives are eliminated. The quantity 6z/60 is zero and 


6z_ = f^a = a o sln ‘V 
6t a dt a sin p 


40 



We can then express the derivatives of c q 


as 


dc 


o 


c)V 


6t u t "St 1 
a o 


a sin 0 
o c 

T 


a cos J 0 


(cos |X sin 7 + sin |x cos 7 cos 0}^r- 

a 


+ (sin (x cos 7 + cos [x sin 7 cos 0} 


d7 

dt 


a J 


+ 


a Q sin7^ ^ 


, du 

da n 

o o 


s in [X \c o s 9 d z 


dz 


28) 


8c -n,, a sin 6 cos [X cos 7 sin 

o _ gV o o 

60 ~ U t Q *3*0* 


cos^d 


(29) 


o 


In equation (28) the derivative of Un with respect to z is 
taken with v constant in the fixed (xp,y]_) coordinate system. 
Thus, du n /dz is to be interpreted as sin v du x /dz + cos v 
duy /dz , °Both du n /dz and da G /dz are, of course, evaluated 
at 8he aircraft altitude only . 

Demonstration of the galilean invariance of A is, of course, 
not essential to the analysis. Here we outline such a demonstra- 
tion. The quantity c is altered by an added constant in a 
galilean transformation in the xp direction, while a constant is 
added to the function u t (z) in a galilean transformation in the 
y x direction. However, the combination of terms 


6t a Ut o dt a 

and the corresponding quantity using 0 derivatives are galilean 
invariant, as are the corresponding derivatives of v alone ^ The 
combinations of integrals I 2 - u t 0 I i and X 3 " 2ut o I 2 + ^o 1 ! 
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are also galilean Invariant. Prom these results it may be shown 
from the expression (26) that A is galilean invariant. 

Application in the program .- For each ray, the derivatives 
of c 0 and v relative to t a and 8 are calculated from ( 27 ), 
(28), and (29). using information from the Aircraft Maneuvers 
Section. The integrals of equations (24) are computed along the 
ray. The ray-tube area A is computed from equation (26) at the 
same time the ray is computed. 

Historical note . - The concept of ray-tube area in sound 
propagation without planar, cylindrical, or spherical symmetry 
appears in a solution by Rayleigh in 1878 (Sect. 284 of ref. 6 ) 
with straight rays. Most discussions in the literature of ray- 
tube area with curved rays have been confined to cases in which 
the aircraft is in steady level flight (no dependence upon t a ); 
in these cases the solutions are much simpler than in the general 
case . 


Flow Near the Aircraft 

The purpose of this section is to define the F-function used 
in the analysis and to present an outline of the local theory near 
the aircraft which leads to the concept of the F-function. The 
F-function is needed as an initial acoustic signal in the basic 
geometric acoustics calculation. The F-function can be directly 
computed by linear theory from the geometry and lift distributions 
of. a slender aircraft. The reader uninterested in the details of 
this computation may skip to equation (33) where the F-function as 
used in this analysis is defined. 

The initial conditions for the calculation of sonic boom pro- 
pagation must be obtained from the flow field near the aircraft. 
This section reviews the local theory near the aircraft which leads 
to the concept of an F-function. This F-function is the function 
in terms of which initial conditions are specified. 

Thus we describe here briefly the linear solution for the flow 
about an aircraft with particular attention to the outer asymp- 
totic form of this solution. We assume, for simplicity, that the 
slender-body and thin-wing assumptions of linearized supersonic 
aerodynamic theory are valid, and that the aircraft may be repres- 
ented by a combination of linear and surface distributions of 
source and lifting elements. 

We assume further, presuming the shape of the aircraft is 
given, that the problem of finding the lift distributions has been 
solved . Thus we shall treat the lift distributions (and corres- 
pondingly the side force distributions) as known. With the slender 
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body assumption, the source distributions are obtainable directly 
from the aircraft shape. For simplicity in notation, we assume a 
single line distribution (fuselage or nacelle) and a single 
surface distribution (the wing). A complete aircraft will require 
several of each type . 

For this section only, we use an (x,y,z) coordinate system 
fixed with respect to the aircraft, with the undisturbed velocity 
V in the direction of the x-axis. The z-axis is vertical; the 
y-axis is lateral. The velocity potential in the undisturbed 
flow is Vx; let V<£ be the perturbation to the velocity poten- 
tial. The reduced perturbation potential O , which has the ^dimen- 
sions of distance, may be divided into the part due to ^ the line 
distribution and the part due to the surface distribution. ^2. ° 
these may be divided into contributions from sources (representing 
cross-sectional area), from lift (z-forces), and from side force 
(y-f orces ) . 

The lire distribution contribution from sources is 


4> 


-1 

2tt 


A ' (x )dx 
v o ' o 

R 


where 


R 


L( x - x o) 


(y - y n ) 2 - " 


O 


\2] i/2 
z o } J 


__ ^2 _ and y = yo > z = Zq is the axis of the line distri 
but ion; the integral is taken to the value of Xq for which R — 0 
to the upstream Mach cone from the point (x,y,z). The lower limit 
-oo simply means far enough upstream to include all disturbances. 
The quantity VA 1 is the linear source strength distribution; 

A(x 0 ) represents the cross-sectional area of the body represented 
by the source distribution, and A 1 (x G ) is its derivative. 

The line distribution contribution from lift is 


3> 


2irpV 2 y 


(z - z Q )(x - x Q ) f z d: 

(y - y Q ) 2 + (z - Zo) 2 R 


and that from side force by the same expression with (z z Q ) 
replaced by (y - y Q ) in the numerator of the integrand _ and with 
f„(x^) replaced by f (x ) . The distribution f^(x 0 ) is the li 

force per unit distance on the axis of the line distribution. 
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Surface distributions of sources and lifting elements are 
distributions on a cylindrical mean surface z = z 0 (x n y 0 ) . T 

contribution to 4> of the surface source distribution " is of the 
form 


t'^o^o) dy Q dx Q 
R 


Here^ t(x 0 ,y 0 ) .represents the wing thickness distribution meas- 
ured in the z -direction, and t' its derivative in the x n 
direction. The lift and side force contributions of the surface 
may oe expressed analogously; here f and f z are replaced by 

distributions of force per unit projected area with an integration 
over y G . 

We next pass to cylindrical coordinates (x,r,0 r ) with 

r = y + z and y = -r sin 0 r , z = -r cos 0 r . With this defi- 
nition, the ray parameter 0 will be given by 0 a + 0„ where 0 a 
is the bank angle. Our purpose is to define for each value of <bt 
a body of revolution equivalent to the aircraft for a distant 
observer. We define further 


r o^ y o’ z o’^r^ = ~ y o sln 0 r " z o cos 0 r 


S o ( ^ y 0’ Z o’ (p r S> = _y 0 cos ^r + z o sin 0 r 


In terms of these variables ^ we can write 

R = [( x - x G ) 2 - P 2 ( r - r Q ) 2 - P 2 s ^] 1/2 


Par from the aircraft, with x and pr both large but with 

x - pr not large, the term s§ will be negligible. We drop this 
term and write 

1/2 1/2 
R = [x - x q + P(r r Q )] [x - x Q - p(r - r Q ) ] 


We next neglect x - Pr - x - Pr 0 in comparison with 8r and 
write 

R = (2pr) 1/2 [x - pr - (x Q - pr Q ) ] 1/2 

In the expression for the contribution due to lift, making the ana 
logous asymptotic approximations, we replace the factor 
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z - z 


) 2 ] 


2 

(z - zo)(x - xo)[(y - y 0 ) _+ ( 

analogous factor for the side force by 


place the dummy variable 


by 


” 1 by -P 
-(3 sin <p r 


cos 0 r and the 
. We also re- 


x 


1 


X 


o 


pr 


o 


In the contributions from the surface distributions the 
integrals with respect to y G are taken with .x-, constan . e 
obtain thereby expressions in terms of new equivalent line dis- 
tributions which are functions of x^_; thus, for example, t 
is defined 




A M (x n ,0 r ) = 


w 


t dy 


x^ =const . 


Distributions f yw and f zw are defined analogously. 

We now define an equivalent area distribution S '(x 1 ,(p r ) 
which replaces all the others, by the relation 


= A' (x, + Pr Q ) + a;( X i ) + 


JL 

P v s 


i(x n ) 


where 


4(x ) = f v (x 1 + Pr Q )cos 0 r - f z (x 1 + Pr Q )sin </> 3 


l A-, 

i 


+ hw 003 h ' sln A 

In these expressions r G is the value for the lane distribution. 
Assembling the terms together gives for 0 the asymptocic ex 
pression 

-1 f Z S ' (X l )d * 1 (30) 


3>(x - Pr J r J ,0 r ) 


2ir(2Pr ) 


172 J (x - Pr - X-, )h 2 


The distribution S'(x n ,0 r ) is the area distribution of an 


equivalent body of revolution on the x axis, as seenby an 
observer a large distance away in the direction . 0^ 
button S 1 (x-, ) is the sum of two terms. The first 


The distri- 
is A' + A^ 


the x 


1 


derivative of the projected cross-sectional area cut by 


the olanes xi — const. .m. — — — 1 - - j.- 

the equivalent force distribution formed by the sum over planes 
x = const. of force components in the direction opposite to < 


The 


second is proportional to i(x 1 ). 


0, 
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With (x,y,z) fixed, the surface R = 0 is an upstream- 
facing Mach cone in the space (x Q , y 0 ,z Q ) . This represents a 
surface of 'coincident signals" for the point (x,y,z). When such 
points are far from the aircraft , these surfaces are approximately 
the planes x 1 = const . (see fig. 9). 

In an actual computation, the determination of the function 
S ( x 1^0r) may be fairly complicated. The force distributions must 
be obtained, of course. Interference contributions generally must 
be computed (ref. 11). If the slender-body approximation is 
inappropriate for any line distribution, additional analysis is 
needed to obtain the appropriate singularity distribution. The 
inlet captured area and the exit jet area must be included with the 
engine nacelle. Finally, the contributions from all line distri- 
butions and all surface distributions are combined. 

With S 1 assumed known, we differentiate equation (30) with 
respect to x . The perturbation pressure Ap is given by 
-pV^<t>x • The result of the differentiation is 



1 C 

2 °p 


_L = ia = 1 p 

M 2 pa 2 M 2 a ( 2(3r ) 1//2 1 


Pr,</> r ) 


(31) 


where 


F.(x 




]_ p S"(x 1 ,0 r ) dx x 

2rr '-I (x - Pr - x-. ) 1//2 

—00 v ^ ]_ 


(32) 


Here we have introduced the magnitude of the perturbation velocity 
q and the perturbation pressure Ap = paq , variables which are 
appropriate for geometric acoustics. Equation (31) gives a solution 
which fits geometric acoustic theory. The r which appears in the 
factor ( 2|3r ) 1 /2 ls a ra y-tube area, x - £r is a phase, and 
and x + P i r ^ are ray parameters. In this case of steady flight, 
the solution is essentially independent of the second ray parameter 
Qtfhich will correspond to t a ). Figure 9 shows the wave front 
x - Pr = const . which are essentially the same family of planes as 
those determined by constant values of the dummy variable X]_ 


Equation (32) defines an F-function in the way 
usually defined. For the purpose of this analysis, 
different definition is used. In place of equation 


in which it is 
a somewhat 
(31) we write 


Ap q 1 

2 _ a ~ 1/2 

pa r 7 


(33) 
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Figure 9. Local coordinate systems. 


The phase variable x - Pr is denoted L , in terms of distance 
along the aircraft axis. This quantity is divided by L a , the 
length of the aircraft, to form the dimensionless phase variable 

k/Lfl . The F-functions are directly related, of course, and we 
write 


P 


Pf P i 


( 34 ) 


The quantity P^ 
given F-function 
With F^ defined 


is simply a conversion factor, to convert any 
Fj_ to the form corresponding to definition (33) 
by equation ( 32 ), this factor is 


P 


f 


2 

M 

(2 p ) 1//2 


( 35 ) 


A plot of this particular function is shown on figure 10. 



Figure 10. F-function factor. 
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between the F-function end weve drag, ^^^“drag represented 

^entrgyTn the^Ivelysrem and ?he rert^drag reprinted by 
energy In tKe trailing vortex system. The wave drag may be 
represented by 2 t t 

(36) 


°D " 


'D 


0 


0i 


d0_ 


where 


C 


D 


0 r 


5 ref 


■ J 


OO 

n 


P, (L,0 r ) dL 


1 


in terras of the F ± defined in equation (32). A transformation 
leads to the alternative expression 


00 x 


C 


8tt^S 


■J J S"(x 0 )S"(x 1 )iin(x 0 - x x ) d Xl 


dx 


(37) 


ref -°° 


Thus there is a direct relation between the F-function and the drag 
of the aircraft . 

The F-function, as we have described it, is a function of the 
ine 1 iunouruii, _ g y an d of the azimuth angle 0 r . 

also e depen!ant upon the^aerodynamic^state of the aircraft^ 
state, wiun some nu^ x ±nclude the dependence upon M a 


maneuver, for example).. We 
in our notation, and write 


M and C 


L 


= F( L/L^ , 0^ > M, Cy ) 


(38) 


are known functions of time t a , as is 


in the remainder of the 
F(L/L a ,t a , 0 ) of phase and 


The quantities M and C L 
the angle of bank 0 a = 0 - 0r • f? us ' 
analysis, F is considered a function 
the two ray parameters. 

The linear theory described in this section does not give the 
only method^ or obtaining ^notions. ^-^ooa^aerod^namro^^ 

the^avfsylteS will yield the P-funetions. They may also be 
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obtained from wind tunnel tests or from flight tests, with perhaps 
™^° r ^ modlflcatlons in the consequent analysis because of nonlinear 

®* am ? le ' det us consider that pressure measurements are 
Wath . a micr °P hone mounted under a balloon in an atmosphere 

W mh5 S ' The aircraft flles by at a distance r in level 
of SKJ pressure perturbation Ap is measured as a function 

is Sen nht^-L/?° m ® a P? earanc< r of the signal. The function F 
hS nf t S° m eduatlon (33). The phase argument of F m 

• f> , is obtained from L = . The nonlinear correction, 

approach Soul”? S aarried ° ut dn two ways. One straightforward 
approach would be to consider the measured F-function to correspond 

init ial C valnp C Slgn 5 ^ has been distorted, to assign an 

fSSthS SSieT°H-°f S® ase ( deflned below in eq. (46)) to account 
sections u?th 1 dj - st ° rtl ° n ’ and fc o proceed as in the next three 

add-on J L T . re ? laced by t - • This approach would entail 

additional analysis to compute x 0 . A simpler method is to apply 

•n e ^ r ° Cei dure of the next three sections to the experimental case 
Tr--p£ Sn?^ Se ' sda rting with the pressure signature and going to an 
i 1 unction . In this reverse procedure, segments of the phase mav 
appear an which P is undefined, or rather is not uniquely 7 
da f a £ ed - nonuniqueness does not matter in applications for 

which the distortion is sufficiently great (t - t > 0) One 
correct way ^ of filling the empty phase segments for the 'purpose of 
his analysis would be by connecting the known portions of the 
curve with straight line segments (with F then continuous). 

Appl ication in the program .- The F-function must be specified 
as anUE-puh function P l( LAa,t a ,0) . The conversion facto? Pr 
is specified or calculated as a function of to . The function F 
analysis° btained ~ rom equation (3^) for use in the subsequent 

tt Historica.1 note .- The theory of this section was given by 
Hayes m 194/ (ref. 12), primarily as the basis of a method of^om- 
putmg wave drag by means of equations (36) and (37). An alterna- 
tive approach was given by . Lomax in 1955 (ref. 13; see also ref. 

Walkden 3n X T > QkR t ^Sf‘ 0f i^? ni S b °? m theory ln a uniform atmosphere by 
/ TJ _ f > d (ref. 15), based on an earlier paper of Waltham 

1 f ' b L l S ° uses an ^-function dependent upon as a para- 
de?!' 1/^ 6 m F - functlon stems a basic pape^by WhithJm 


Geometric Acoustics and Blokhintsev Invariance 


The purpose 
iant measure of 
valid globally. 


of this section is to define an appropriate invar- 
the. acoustic signal along each ray, one that is 
This measure of the signal will be expressed in 
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terms of the F-function for the aircraft. 

reometric Acoustics is a linear acoustic theory based upon 

the assumption that _ the solution a )f The wave 

si'D‘D6cirs locally as i*t would in a plan If , ' -.i . -Hopv 

???n?rare actually gently curved n -^ es ih ^* f ^ C e f gkea ?he 

pressure Ap given hy 

Ap = paq (39) 

As in a plane wave system, the quantity q is a function of phase. 
When the undisturbed flow field is steady in a specific ^ 

Ssl?vat^ me ^f f S°phasS ^The^efer- 

found in reference 4. 

The fact that a birectly a defined f , r useful, ^sical^ntit^ ^ 

restrictior^tif steady atmospheric properties^^^This aa n ^anrent ^ 

element would not appear in an acous t hyslcal inter- 

atmosphere, where a phase variable with 
pretation would have to be used. 

The ratio q/a is. used as 

In a plane-wave system in a u nifo g ach wave front and thus 

he a function only of With the atmosphere non-uniform 

constant on each kinematic r y. no longer constant on 

and with th ®. wa ^® fr0 5^ S 0 ne r Ipproach to this intensity problem, 
each kinematic ray. In PP • j- eg . ra ]_ along the ray (refs, lo 

the intensity is expressed a J ^ r g 1 | includes a term 

and 19 ). m this approach the ^tegrand ax y must then be 

ISSSS^^oSS-^tSSSl^e alternative approach 
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measures d ®P enda u P° n theoretically established invariant 

rty-tnlt area intensit y ■ This approach requires the concept of 

the Revi Scb t a? SPh ^ S at rest ' wlfch n ° winds, volume integrals of 
tne Rayleigh acoustic energy density (ref. 6, Sects 24s and 

are conserved. _ This acoustic energy density'is p^/1^%^ 

Within geometric acoustics, with equation (39) valid, this enerev* 

potential 3 B ^ Ply pQ 1 Half the ene ^ ^ kinetic hSd hal? 87 
1 * ^ he en f r Sy flow down a ray-tube is pq 2 aA^ , where A 

ihus i?° r n/e ray l Ube area ^ md is constant on a kinematic ray > 
inus ii q/ a is known as a function of phase at one point on a 
geometric ray, it can be predicted at another. Rayleigh (ref 6 
e?S; 266 and 284) used essentially this approach to predict ‘the 
effect of density and ray-tube area changes. 

In an atmosphere with steady winds, an analogous invariance 

?ef PIT 8 f? U 2 d h B1 ? Ehlntsev ln 1946 (Sect. 7 of ref. 20, or 
r;? 1 * I • abld only Within geometric acoustics. A straightfor- 
w a rd derivation may be found in reference 22, Sect 1 The invar- 

^ the quantity that is'cltant a?Sng 
Kinematic ray is this density times the volume flow c . A Me 

V S C °k^ a ?i quantit ? as e D times the square of a ' 
hevieh J E ^ 4he_sign of V E the same as that of q We 

have already defined A as A times the unit vector in the -2 
direction. The invariant is thus 


C 0 V E 


22 - 


Ac a sin 0 

ct n 


or, using the relation 


n 


o 


cos 0 (eq. 14), 


V E ( *1, t a ,$) = ( pa 2 A sin 0 cos 0) ^ 


(4o) 


c 


factor in the definition of vi 


E 


is to 


The purpose of the 

^ aE ® bhe s Y bs ®4uent analysis simpler and also to make ~Vtu Gali- 
lean invariant. Once V E (|) is known for a particular geometric 
y, we can solve directly for the pressure perturbation Ap(£) as 
/h n V nc aon °T the actual time t + % using equations (39) and 

tion * fo^fhp t- 1S ^ tlme obtalned ln the ra y tracing computa- 
tion fop the kinematic nays of zero phase . 

General DejtinH^ n? tC \ th ? crltlca l step mentioned in the 

a of the analysis, in which a consideration of 

a galilean transformation is inescapable. The local solution is 
given as a function F of L/L a and two ray parameters in a 
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coordinate system moving with the velocity of the aircraft. The 
global solution is of the form of a function V E of £ _ and the 
two ray parameters t a and 0 in a coordinate system fixed rela- 
tive to the ground. In order to express the solution in terms of 
V e (4) > we must relate both the dependent and independent variables 

in the two coordinate systems. Although what has technically to 
be accomplished is the galilean transformation, much of the job has 
already been done by the identification of ray parameters and phase. 

The ray parameters give us no trouble, as the galilean trans- 
formation has already been accomplished in the details of calcula- 
ting the initial wave normals and of the ray tracing. .The 
parameter t a is identified locally through the functions M(tg,), 

C (t ) , and 0 a ( t a ) • The local and global azimuth angles are 

connected through the relation 0 = 0a + $r • 

To relate the two phase variables, we use the stratagem. of 
introducing a third, local, phase variable which is readily inter- 
pretable in either the local or ground-fixed coordinate system and 
has the property that it is galilean invariant. Such a variable is 
the distance s normal to the wave fronts, measured from the 
reference wave front at a given instant. In the local coordinate 
system, we have 

s = L/M = L sin p 

The quantity L is distance from the reference wave front in a 
direction \ tt - p from the normal. In terms of the phase £ 
introduced in this section, the distance s is 

s = ^c no = £o 0 cos 0 Q 


We equate the two expressions for s and obtain the basic phase 
relation 


L a sln ^ L 

c cos Q~~ L 
o o a 


( 41 ) 


The presence of the factor c Q in this relation indicates that £ 
is not a galilean invariant quantity. 

To relate the dependent variables V E and F , we follow an. 
analogous course. The dimensionless measure Q./a of the signal is 
clearly galilean invariant and may be expressed in terms of either 
F or V F . We equate the two expressions for q/a . This relates 
the dependent variables, of course, but the connection is incom- 
plete. The ray-tube area is singular at the aircraft, is repres- 
ented by the variable r in the local coordinates and by A in 
the ground-fixed coordinates. We must then also relate the two 
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corresponding ray-tube area measures to a measure which is clearly 
galilean invariant. Such a galilean invariant measure is time 
(t - t G ) elapsed from emission of a signal from the aircraft axis. 
In the local coordinates of the previous section, the quantity r 
is related to this time by 

r = a Q cos p (t - t Q ) 

The definition (33) of the function F leads to 


q 

a 


[ a Q cos p (t - t Q )] 


1/2 


(42) 


For the function A we must carry out a local calculation of 
equation (26) near the aircraft. To lowest order we have 


A = if 1 + — r- sin 7 \ cos 0 Q V sin 0 cos p cos 7 

L\ sin M- sln e 0 /sin ~ q ^0 + " 


cos 0 


Q 

COS J 6 


X 


o 


a sin^0 
0 o 


/ dc 
l 0 

6v 4 " 

\ 60 

- U f ) 

to 60 / 


(z - z ) 
v o ' 


The quantity (z Q - z) is approximately equal to a sin 0 o (t - tp). 
The other terms may be re -expressed with the aid of equations (10), 
(27), and (29). After some calculation, we obtain 


A = 


2 

a cos p cos 0 (t - t ) 
o o v o ' 


sin p sin 0 


o 


From equation (40) we obtain 


q 

a 


, x 1/2 

(sm p) V E 


[p a-( cos^p cos'^0 (t - t )1 
L *0 o ^ o v o ' J 


172 


(^3) 


( 44 ) 


Equating the two expressions (42) and (44) for q/a gives 
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F 


(45) 


V. 


2 „ 2 fl 1/2 
,/p a_ cos M- cos y , N 


E 


o o 


o 


sin. (i 


This Is the desired relation, giving V E In terms of the F- 
functlon. 


Application In the program.,- For each ray Investigated, the 
function V F is computed from equation (45). This is expressed 
in terms of the independent variable £ obtained from equation 

(41). 


Historical note.- The invariance result of Blokhintsev was 
found independently by Chernov in 1946 (ref. 23), but only in the 
special case of ipz^ofcational flow. Gar*]?ett (r*ef . 2m-) noted, that 
Blokhintsev ' s result was equivalent to conservation of volume 
integrals' of pq 2 /f2 where pq 2 is the Rayleigh acoustic energy 
and H a frequency, both measured by an observer moving with the 
fluid. Hayes (r ef - 2 5> 1968) showed that Garrett's result remains 
valid when the undisturbed flow_field is unsteady, and that the 
quantity pq^aA n /2^ = pq^c n c • A/o> a is constant on a kinematic 
ray in this case, where A n is a ray-tube area cut by wave fronts 
and cd is a frequency measured by a fixed observer. Bretherton 
and Garrett (ref. 26, 1968 ) present a general theory governing wave 
motion in moving media. 


Signal Distortion and Age Variable 

The purpose of this section is to define an "age" variable and 
to apply it in the calculation of the nonlinear distortion of the 
acoustic signal. This distortion appears in the propagation of 
acoustic signals over large distances. The distortion's caused by 
a weak, nonlinear effect resulting from small changes in propaga- 
tion speed which are proportional to the strength of the signal. 
Although the nonlinear effect is locally weak, it is cumulative and, 
as a consequence, the total distortion of a sonic boom signal is 
far from negligibly small. Shock waves may appear in the signal 
or two shocks (or more) may merge into one. The process of distor- 
tion is governed by an age variable which is defined and applied in 
this section. The study of the location and motion of shock waves 
is to be covered in the following section. 

We make here a simplifying assumption, one made by all invest- 
igators of sonic boom. We assume that in a lowest-order approxi- 
mation the phase shift due to the change in propagation speed is 
the only nonlinear effect that needs to be considered. Thus, any 
nonlinear effect on the rays, the ray-tube areas, or the Blokhint- 
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sev Invariance is assumed negligible. Although no comprehensive 
theory is available to fully justify this assumption, it is not 
made blindly. The neglected effects may be shown to correspond to 
higher-order terms in a small parameter for level flight in a 
uniform atmosphere (re'f. 27 , section on Finite Systems). In 
general, it appears likely that this assumption fails to be valid 
_only when the assumptions of geometric acoustics fail or when the 
perturbations are no longer weak. Nonlinear effects on rays are 
discussed by Whitham in reference 16 and form an essential part of 
his method for predicting the trajectories of finite strength 
shocks (refs. 28 and 29). 

The propagation velocity, equal to c = an + u in the undis- 
turbed fluid, is changed by (Aa + q)n , where Aa is the pertur- 
bation in the speed of sound. This quantity is given by 


Aa 


7, 


- a 

2 pa “ 2 q 

so that the change in propagation velocity is simply (7 + l)qn/2 . 

If the phase were expressed as distance s measured back- 
wards normal to the wave fronts, the signal would experience a 
phase shift arising from the change in propagation velocity given 
by 

ds 

dt ~ 


7 e - 1 


7 e + 1 


q 


where t is time along the ray, that given by equation (l8). 


In treating the phase, we must distinguish between the actual 
phase variable and the phase variable according to linear theory. 
The nonlinear effect is basically the difference between the two. 
We term the actual phase ^ , defined in the same way as before, 
and term the linear phase 4 • The local distance phase s is 
related with |-j_ by ds = c n d|-|_ . The expression for the change 
In propagation velocity in terms of |-j_ becomes 


d| x 

dt 



1 




1 


q 

C Q cos 9 


This equation describes the phase shift for a particular point on 
the signal, and thus with the linear phase | fixed. The point 
on the signal found at | according to the linear theory is 
actually found at i 1 , which is a function of | and t . 

We express q in terms of using equation (40), and 

obtain 
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y +1 

'e 


cos 0 ( pA sin 0 cos 0) 


The variable t may be replaced by -z through equation (l8). We 
wish to transform the equation for the phase shift to a canonical 
form. We introduce the age variable r defined by 


y +1 

r e 


d( -z 


a sin 9 cos 9 ( pA sin 9 cos 9 ) 


(46) 


The phase shift is then governed by 


V E (q,t) 


in the canonical form desired. 

With our basic assumption that the linear results for the 
ray-tube area and Blokhintsev invariance still hold with only 
the phase shifted, we must have 

V E (e i’ x) = V E (I) ( 

where Vv ( £,) is the linear solution (independent of t or t 
because of the invariance). Here the actual phase S4 
satisfies equation (47) at constant | and also the inxtia 

*• ~ The solution of equation (47) 


condition 
is then 


— 4 at ^ 


= £ ~ ^( 1 ) 


(49) 


The linear phase i is a function of ^ and t which may- 
become multivalued in . In this case the solution ( ) 

and (49) must be modified to take into account the presence ox 
shock waves. This modification is explained in the following 
section . 

The solution (48) and (49) is the . general solution to the 
first-order partial differential equation 


(50) 
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This equation, may be considered a canonical one Tor inviscid wave 
propagation in one direction. One approach to the problem of 
calculating the nonlinear distortion lies in first deriving an 
equation of the form of equation (50). Equation (47) or its 
integral (49) gives the characteristics for equation ( 50 ) and its 
solution by standard methods is that of equation (48). 

If V E is set equal to the 4q derivative of a function xf 
in equation ( 50 ) and the resulting equation is integrated with 
respect to £ >1 , the result is 


cb/_ _1 f 

6 x 2 \c)^/ 


(51) 


We have here assumed that the function is zero for sufficiently 

large negative with the consequence that the additive arbi- 

trary function of t from the integration is identically zero. In 
the following section, we define the function as an integral 

over | and use its properties to locate shocks. Equation ( 51 ) 
may be considered as an equivalent of equation (50). Actually, 
because of its properties when there are shocks, the function 
contains more information than does V„ . 

ill 

If the P-functions are obtained experimentally from wind- 
tunnel or flight test results at distances from the aircraft at 
which x would not be negligible, a correction may be needed. 

This could involve a correction transforming F to an equivalent 
F as mentioned in the section. Flow Near the Aircraft. An alter- 
native method would require calculating an initial value of x in 
equation (46) which is not equal to zero. 

If t gets very large, the signal shape approaches that of 
an N-wave. With the shocks properly taken into account, the 
asymptotic behavior includes a maximum to |V E | which is propor- 
tional to t / and a total phase difference between head and 
tail shocks which is proportional to x 1 / 2 . The presence of the 
factor p V 2 in the integral of equation (46) indicates that for 
downward propagating rays in a real atmosphere the integral will 
be convergent; thus x approaches a limiting value x-, • as z 
approaches -00 . The signal in terms of V E approaches the 
limiting shape VE(4]_,X]_i. m ) . Thus the signal shape "freezes" 
and does not become the ever-thickening N-wave of common farfield 
theory. This freezing effect is discussed in reference 1 . 

To illustrate the freezing effect a simple example is 
instructive. With level flight in a constant temperature atmos- 
phere the asymptotic value of the age in a downward 6 > Iv) 

direction is proportional to 
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00 / p \l/2 

J ( z o - z)" l/2 exp[ -|p(z Q - z)] d(-z) = \jr) 


-z 


wher 


-*e p" 1 is the scale height of the atmosphere. 

The corresponding finite integral in a homogeneous atmosphere 


is 


(z - z ) 
v o ' 


- 1 -/ 2 d(-z) = 2(z o - z) 1/2 


-z 


o 


and the two integrals are equal if z Q - z = m/2p * J; h ® 
signal shape in the real atmosphere is, therefore, the same as 1 
would be in a homogeneous atmosphere at an altitude 7r/2 scale 
heights below the aircraft. 

Applicat ion in the program .- For each ray, the variable t 
is calculated from equation (4oj at the same time the ray and r y 
tube areas are calculated. At the ground the original Phase £ is 
transformed to the actual phase ^ by equation (49), giving 
SS?eby the function VeO^t) at the ground. No provision is 
made for the correction required if the F-function is obtaine 
from a measurement made a large distance from the aircraft flight 

axis . 

Historical note.- The solution (48) and (49) was obtained by 
Poisson in l80Y (ref. 30) for plane waves in a constant -temperature 
gas (7 =1). He obtained also equation ( 51 ) corresponding to this 

solution , with J a velocity potential ( ^hich 

31 ) showed that with a polytropic gas the factor 2 (7 e + 1) wnic 
appears in equation (46) must be included m the analysis. An 
equation of ?he form of equation (50), although equivalent to 
Poisson's integral (51), seems not to have appeared before the 
issuance of a report of Chandrasekhar m 1943 (rel . 32 ). 

With planar waves in a uniform gas, the age variable is simpiy 
proportional to the distance coordinate. The . simplest cases wit 
variable ray-tube area are those with cylindrical, spherical, or 
coSXl symmetry. In these oases, for N-waves in a uniform atmos- 
phere , Landau (ref. 33) obtained the correct laws Th solution 
with a general signal, corresponding to that of this and the 
following section, was obtained by Landau and Lifshi irh^-f-ham 

(ref. 34) with cylindrical and spherical symmetry and by Whitham 
in 1952 (ref. 17) with axial conical symmetry (flow about a body 
of revolution). The application to flows with general conical 
symmel?y? with the azimuth 0 as an independent parameter, was 
Sp hv Haves in 1954 (ref. 27, section on Finite Systems). The 
arbitrary dependence of the solution of 4 here involves the 
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principle that nonlinear effects on the rays are negligible. Ray- 
tube areas not based upon solution symmetry appeared in the age 
variable used by Rao (ref. 35, 1956) for maneuvering flight in a 
a ^ raos Pfr ere with straight line rays. Here the assumption 
that nonlinear effects on rays are negligible is implicit. 


Variable atmospheric properties appear in the age variable 6 
used by Whit ham (ref. 3 6, 1953) with spherical symmetry in a study 
of weak shocks in stars. Arbitrary area and fluid properties were 
combined. by Hayes in 1957 (ref. 37) without consideration of winds, 
several incorrect age definitions appeared in the Russian litera- 
ture for the case with steady winds. In 1962 Ryzhov and Shefter 
/£■ us ® d , the correct one but applied only to N-waves. Hayes 

’ d 9o3j presented an analysis with winds similar to that of 
this section but made an error in the definition of the age (a 
factor a/c n is required under the integral); this paper includes 
an algorithm for ray-tube area with winds. Guiraud (ref. 19 1965) 
uses an age variable* but one defined somewhat differently because 
he does not use Blokhintsev invariance or* nay-tube areas , 


The freezing effect , although inherent in the age variables 
appropriate with varying atmospheric density, appears not to have 
been discussed before reference 1. The integral involved with 
level . flight in a constant temperature atmosphere is an error 
function and appears in a related form in the appendix to a 1955 
paper of Busemann (ref. 39). 


Shock Location 

Tne purpose of this section is to show how shock waves which 
appear in the distorted signal may be located. In the algorithm 
this is done at the same time the distorted signal is calculated. 

In general, with the age t sufficiently large, the distorted 
signal V E (£]/!■) is multivalued in ^ and, thus, physically 
meaningless. The actual signal contains one or more shock waves, 
and a proper treatment of the shocks eliminates the multivaluedness. 
A shock wave moves faster than the acoustic propagation speed in 
front of it and slower than the propagation speed behind it. 

Thereby, parts of the signal are propagated into the shocks (are 
eaten up by the shocks) and this phenomenon permits the remaining 
part of the signal to be single -valued . 

The procedure for locating the shocks utilizes the function 
<*7 (Ii,t) . introduced in the preceding section. We drop the t in 
the . notation and write simply >/ ( ^ ) . in order to have 
defined over the entire range of 4 > it is defined as an integral 
over | , or, equivalently, as a Stieltjes integral over £, (£) 

Thus we write . 
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( 52 ) 


= i(«) -/var 1 )^ 

— OO —oo 

This function satisfies equation (51) over the entire range of | . 
The lower limit -oo simply means sufficiently negative to begin 
before the signal commences. The portions of the original signal 
Vg(|) which appear in the actual signal correspond to segments or 
branches of the multivalued function >f(|q) , chosen in such a way 

as to make up a single-valued function. 


The bisector law for a weak shock states that it moves with a 
normal speed midway between the acoustic propagation speed in front 
of and behind the shock. This law may readily be proved, and 
requires only that the curvature of the fluid isentrope be contin- 
uous. In terms of our notation., the law states that 



shock 


- 5< v E 


+ V E ) 


+ 


(53) 


where + and - indicate the two sides of the shock. We use [ ] 
to denote jump in a quantity; thus [>ft = < -yf. ■ The t 
derivative along the shock of [*f] is calculated using equations 
( 51 ) and (53); and is 


dT 


~cb/ 


_dx_ 



.dr j 


©. 


shock 


1 


— ( 

2 ^ E 


- o + 


d|- 

dx 


■) 


shock 


(V E - V E ) = 0 


E 


The function -J is continuous in the original signal at x = 0 , 

so that [?f] is initially zero at the point of formation of any 
shock. Thus, 3/ remains continuous across any shock. .It is clear 
that this property is preserved when two shocks merge into one, 
and that the property thus holds in general. 


We conclude then that the single-valued function of £q made 
up of segments of |q) is a continuous function. This result 
we may term the equal-area law. In its simplest form the law .states 
that a shock, which appears when the final curve. jumps from one 
segment of -d ((=i) to another, is located where it cuts off equal 
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Two important properties of >/(! i) are not evident from the 
analysis above, and we state them without proof. One property is 
that there is only one continuous function made out of segments 
of «/(£-,) for which d^/dS >0 . The principal consequence of 

this property is that we do not have to trace the trajectories of 
the shocks from t = 0 but may locate them at a fixed value of t 
without considering their history. The second property is that the 
desired function is precisely sup*/(!q) , the superior limit of 
/(g ) defined to be the function obtained by taking the largest 
value of J' available for each value of ^ • This property sim 
plifies significantly the process of locating the shocks m a com- 
plicated case. 

The equal-area law may be established in other ways than that 
of using the bisector law. One alternative approach invokes the_ 
principle of conservation of mass, identifying yJ with a Lag rang i an 
variable or with a particle displacement variable. Another approach 
involves establishing equations with viscosity and using the limit 
III violets V i 0 . In the second of these alternative approaches, 

the property that the desired function is sup^d^ appears 
naturally . 

The process of locating the shocks may now be described. The 
function >/(£-,) is obtained, and the single-valued function 
sup sf(Z -i) is 1 identified. The function V e (£ 1 ,t) is then plotted, 
retaining only those segments of which appear m. sup;/. The 

corresponding function Ap(|t,t) is calculated and gives the esire 
pre ssuresignature . The parts of the original signal corresponding 
to segments of not represented in sup >/ are the parts eaten 

up by S the shocks. 1 Where there is a jump _f rom _one segment o£ 4i to 


another in sup >/ , there is a corresponding jump in 
and such a jump represents the shock. 


V- 


E 


and Ap 


The function Wdl) may be obtained by direct integration, 
but a somewhat different construction is a « v ?htageous. Together 
with the original function V E (£) , we require its integral 


s o U) 



(54) 


From equation (49) we have 

6i dV E (l) 

d| 1 ~ T “dl 

at constant t . Applying these relations to the integral (52), 
we obtain 
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(55) 


v^i) = S 0 (0 - I T V E (?) 2 

with equation (49) giving the relation between and f 

Equations (49) and (55) may be combined to give the relation 

■s/^i) = S q (£) — (56) 

The graphical method of Burgers (ref. 40) is a method of finding 
sup ,/(£]_) from a given "summation" curve S Q , using equation ( 56 ). 

The pressure signature Ap is to be calculated from V R using 
equation (40) and the relation Ap = paq (equation (39)). We are 
usually interested in the pressure at the ground, and must take 
into account . the reflection of the wave system from the ground. 

This reflection may be visualized as coming from the mirror image 

with respect to the ground of the' impinging wave system, and the 

reflected pressure signature appears superimposed on the incoming 
one. Very near the ground, say at a height h , the reflected 

pressure is the same as the incident pressure, with a phase delay 

(in ?]_) equal to 2h/a sin 9 . If we assume that h is small 

enough that . the phase delay may be neglected, the effect of the 
ground is simply to double the pressure. This factor appears as a 
reflection factor equal to 2 

In practice, an empirical reflection factor K R other than 
2 is often used. The final expression for the pressure signature 
from equation (40), is then 

Ap(l 1 ) = K r V E (l 1 ^r)(pa 2 /A sin 9 cos 9) 1 / 2 ( 57 ) 

The pressure signature far from the ground, with no reflection taken 
into account, is obtained from equation ( 57 ) simply by setting the 
factor K r equal to 1 . Our analysis and algorithm for cal- 

culating sonic boom signatures are completed by equation (57). The 
ray tracing and ground intersection calculations discussed earlier 
tell us where and when the pressure signatures occur. 

The pressure on the ground at any instant is approximately 
constant along the ground intersections of the wave fronts ( the 
y -j_ direction), and has its main variations in the direction normal 
to these ground intersections (the x-, direction). Distance in a 
plane z = c ons t . normal to the wave front intersections (in the 
X 1 direction) from the zero phase wave front may be used as a phase 
for describing the pressure on that plane at a given instant. This 
distance phase variable is equal to c 0 ^ . This phase is galilean 

xnvariant; it was not used in the analysis as the basic variable 
because appeared to be a more useful quantity. 
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Application in the program .- For each ray, the function 3n(4) 
is obtained from equation (p4)7 and the multivalued function >A?l ) 
f* r>nm prmations (55) and (49). The branches corresponding to 
sSvtiP are identified, and ApUq) is calculated from equation 

(57). 

Historical note.- Shock waves remained ill understood through 
almo st lir - ST^ nTneteenth century. The bisector law was given 
bv Crussard in 1913 (ref. 4l), who described what we would term a 
hal?T£ave In a constanfc-ara duct. This law was used ty Landau 
and Lif shitz (ref. 34, 1944) to obtain the equal-area law (see also 
ref 17). In 1950, E. Hopf ref. 42) gave a thorough study of 
Burgers equation with particular attention to the limiting process 
u, — *■ o . This study establishes the equal-area law through this 
limiting process and serves as the basis for Burgers graphica 
method (ref. 40). Landau (ref. 43., 1945) mentioned conservation of 
mass as the basis for the equal-area law, but without giving 
details. Ligh thill (ref. 44, 1956 ) noted that the equal-area law 
is equivalent to continuity of a Lagrangian variable. Middleton 
and Carlson (ref. 45, 1965) applied the equal-area law m prac- 
tica^calcuiations, obtaining^ j/(q ) through direct integrations. 


A Note on Viscous Effects 


Although they are not properly part of this . analysis , this 
exposition would be incomplete without some mention of viscous 
ISecis In sonic boom problems of practical interest viscous 
effects ‘are sufficiently weak so that their only effect is to 
give finite thickness to shock waves. This . thickness, thoug. 
finite is generally several orders of magnitude smaller than any 
of the* other characteristic scales of the problem. Hence the 
treatment of shocks as strict discontinuities is completely soun . 
Tn Stain other problems, as for example of acoustic propagation 
a? ionosphericheights , viscous effects are important. Here we 
indicate briefly the governing equation, and describe a few of 
its properties. 

As in much of our development, we take advantage here of the 
fact that the wave system is approximately a plane wave system. 
The viscous effects are considered to comprise only viscosi y 
and heat conduction. They appear, in a coordinate system (s,t) 
moving with the undisturbed inviscid propagation speed, m a. 
diffusion term of the form 

^ V E _ _i/4 ^ + p,« + (7 e ~ (58) 

" ^pV3 ^ ^ c p hs 2 
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1S the distance coordinate which serves as the phase 

and ^ are the shear and dilatational viscosity 
• ° e f5 1Cient ?j,. is fche ratl ° of specific heats of the gas, c 

is the specific heat capacity at constant pressure, and k is tlhe 
heat conduction coefficient. The variable P V E in' elation (58) 


^ . x.uc vcti-ictoie vtt m eauation ^ 

may be any intensity measure, but it is identified here as the 
defined earlier. The nonlinear term corresponding to that in 
equation (50) has been omitted. 


V- 


E 


We next change the variables t and s to 

variables x and j»i and include the nonlinear 
result is 


the corresponding 
term. The 


where 


dV E dV K 1 

■3V = V E 3 ^ + I N(t) 




N 



(59) 


is a. positive reduced kinematic viscosity 
is simply c n = c Q cos 9 , while 


The quantity ds/d^ 


dt _ 2c o cos 9 (P A sln e cos 9) 1 / 2 
dx I7 e + lj 


from equation (46). We obtain thereby 


N = 


_ 2(pA sin 9 cos 6 ) ^Z 2 


(7 e + 1 )P C 0 cos 


P + P 


( 7 e - i)k' 


(60) 


that^he auantltv a f T ene J al j zed Burgers equation, generalized in 
onar toe quantity N may be a general function of x 

Lighthill (ref . 44) introduced the concept of the Reynolds 
number of a lobe of a solution of Burgers' equation As aoDl^ed 
to equation (59), we identify a lobe Is that^rt 4 the signal 
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between two consecutive zeros of Vg , ^1 ^ a -, nhps'^can 

Lobes can disappear with increasing t , but no new lobes can 

appear. The first moment % 1 of the lobe is 


b 


foA-c) = 


V E 


a 


The Reynolds number of the lobe is 


Re(t) - 


2 7 % 


1 


N 


(61) 


If the Reynolds numbers of all the main lobes of a signal are very 
larve viscous effects are unimportant and only determine . shoe 
thicknesses. If the Reynolds number of the largest lobe is smal , 
nonlinear • effects are unimportant. 

The first moment of a lobe obeys the equation 


--i = •asr 

dr 2 


dv. 


E 


3? 


r b 


Sv e\ 


C>2 


i* .r sy* ;r,ss „!?ii as rst-ary 

which V E changes sign. The Reynolds number can increase, but 
only if N decreases sufficiently fast with t 


The second moment of a lobe is defined 

>b 

’1 


^2 =jf V E d? : 


and represents the wave energy in the lobe. It satisfies the 
equation 


C 

dr 


= - N 


pb /c>V \ 

L GrJ 


(63) 
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This moment cannot increase. It is constant for 
there are no shocks in the lobe. These results 
paper for which reference 37 is the abstract. 


N = 0 only if 
were given in the 


Application in the program . - 
in this program. 


Viscous effects are not included 


Historical note - The combination of coefficients which 
appears m equation (58) and which governs viscous weak wave propa- 
gation appeared m a 1910 paper of G.I. Taylor (ref 46) on the P 
° f Weak shock waves (with |i' =0). Lighthill (ref 44 
1956) showed that Burgers' equation was the appropriate one for* 
general Planar weak wave propagation for a calorically perfect gas 

Ha|^ e fref 2e ^7 B 4iV! e VV4° n ^ the f °™ ogtSnedby® 

nayes (ref. 37, 1957) without winds. Hayes (ref 47 iq^ft) si^o 
showed that the combination of coefficients in equation ( 58 ) was 

SK: W VL a . g T ral ^r t±0n ° f state « for a WaS 

?Iq\ * S 5 - ? he inclu sion of winds in a derivation of equation 

froi ne by 5 ay ? S (ref - 38 ' 1963 > but w3th error coming 

from the error made in the definition of t An equation of fW 

eqSivileVo? ?h a 4 ed by Gu , iraud («*• b9 , i9 65) inttrZ % ITs 
of visoouS variable V E . . Further details and discussion 

viscous effects may be found in references 44 and 47 


bummarizing Statement 

-p ^ At s poln V the analysis is complete and ready to be used 
for computing sonic boom. Here we briefly summarize the algorithm 
for obtaining pressure signatures in the form realized in the 
computer program. 

atmospheres 1 " 6 S ° nlC pressure signatures in a stratified 


1. Specify input data. 


a. The _ thermodynamic properties of the atmosphere and the 
horizontal wind velocities are needed as functions of the 
altitude z and are used throughout the calculation. 

b. The aircraft Mach number M , its heading angle f , its 

angle J > lts bank angle 0 a are needed as 

unctions ^ of the time t a . The initial location of the 
lr craft is needed to start the trajectory calculation. 

As an option, aircraft load factors may be specified as 
functions of time t 

cl 


c. F functions for the aircraft are needed in sufficient 
number to serve over the range of aircraft Mach number 
covered in the trajectory. 
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2. Calculate aircraft maneuver functions using lb. 

a. The aircraft trajectory, x a (t a ) , ya(ta)> and z a( t a)^ is 
calculated . 

b. The derivatives of heading, climb, and Mach angles, dt^/dt a , 
d 7 /dt a , and dp/dt a , are calculated, either directly from 
the input data or (optionally) from aircraft load factors. 
These derivatives are needed in the ray-tube area calcu- 
lation 4b. 

3. Calculate the initial wave normals using lb. 

The orientations of the initial wave normals, parametrized 

by their azimuth angle 0 , are needed to determine the 

invariants c Q and v used in the ray calculations 4. 

4. Calculate rays and functions along rays using 3. 

a. The rays, described by functions x(t a ,0,z), y(t a ,0,z), 
and t(t a ,0,z), determine where and when the sonic boom 
signal hits the ground. 

b. The ray-tube area A is calculated along each ray using 
2b. This ray-tube area is used in calculating both the 
age 4c and the final pressure signature 6c. 

c. The age t is calculated along each ray. This is used 
in calculating the signal distortion 6. 

5. Calculate the linear acoustic signal on each ray. 

The function V E (4) is obtained from an F-function by a 

transformation of both dependent and independent variables. 

The integral S 0 (|) of V E is also calculated. 

6. Calculate the distorted pressure signature. 

a. For the age at the ground, calculate the distorted signal 
V e (4]_,t) and its integral xf . 

b. Locate the shock waves using >f . 

c. Calculate the final pressure signature Ap(^) using a 
selected reflection factor. This is the desired output. 

The logical order of the algorithm follows closely that of the 
theoretical analysis, with one exception. The age is computed at 
the same times as are the rays and ray-tube area rather than later 
with the distorted signal. 
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COMPUTER PROGRAM 


General Description 

The computer program for solving the preceding sonic boom 
equations has been written in FORTRAN with the aim of obtaining a 
versatile computational technique adaptable to a variety of 
computers. This program was developed on an IBM-1130 and has also 
been run on a CDC-6600. A complete program listing for the CDC- 
6600 is given in Appendix A , and a sample printout appears in 
Appendix B. This printout will be discussed in detail in the next 
chapter, COMPUTATION RESULTS. 

The organization of the program is shown in figure 12. The 
input processor reads input data and converts them to appropriate 
units and format for use in the subroutines. The input data are 
then listed on the output sheets (Appendix B) for identification 
and checking purposes. The program next computes and lists out 
pertinent maneuver data at the selected initial time t a . These 
outputs include aircraft location, Mach number, direction angles, 
and derivatives. The program then proceeds to the ray path and 
ray-tube area calculations beginning with the ray parameters 0=0 
and t a . The program computes the ray trajectory, angle with the 

horizontal (cos 0) of the wave normal, the ray -tube area, and the 
age variable. These outputs are listed as functions of altitude 
z ; the time and (x,y) coordinates are cumulative from the first 
maneuver point. When a stopping condition for the ray trajectory 
is reached, such as z = 0 , the program computes phase variables 
and the -function, and utilizes the F-function to determine the 
distorted pressure signature. These output quantities are then 
listed in the order of increasing values of L/LA . At each ray 
intersection with the ground, the program stores parameters which 
give the location, time, and maximum pressure of the calculated 
signal. The program then returns to the same maneuver point (t a ), 
increments the ray azimuth angle 0 , and proceeds again through 
the ray tracing and pressure calculations. When calculations for 
the set of rays are completed corresponding to the maneuver point 
t a , the program proceeds to the next maneuver point t a+1 and 
the ray tracings are repeated. In the event the ray-tube area 
diminishes to zero during its traverse, or if the ray becomes 
horizontal, the calculation along that ray stops. After completing 
all of the maneuver, ray tracing, and pressure sequences, the 
summary data which have been saved for each ray are used. Linear 
interpolations are made to determine the locations of the rays and 
maximum pressures at the same elapsed time t . Thus, ground shock 
intersection coordinates are made available, together with indica- 
tors of boom intensity there. 
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Subroutines and Operations 

The program consists of a master program, SONIC , and. nine- 
teen subroutines. Figure 13 gives a representation of the re- 
lationships among these subroutines. A brief description of 
their functions is given in the following paragraphs. Numbers 
correspond to those shown on the flow diagram. 

1. SONIC This is the master program which initializes and 
calls the subroutines. It also serves as the main program for 
the maneuver calculation. 

2. INPU This subroutine reads in all input data, checks them 

for self-consistency and validity, and prints out a copy of the 
data for comparison and. problem identification. 

3. ATMO This program calculates atmospheric data when using 

the standard atmosphere option. For a given altitude z , the 
corresponding temperature, pressure, density, and speed of sound 
are found. The altitude above sea level must be less than 

170 600 ft (52 km) . 

4. ANGLE This routine corrects input data so that no two con- 
secutive angles are more than l 80 ° apart; for example, 25 ° and 
- 190 o would be changed, to 25 ° and 170 °, respectively. 

5. DERM2 This routine calculates the integrands used, in 
solving for the maneuver path. 

6 . 0UTM2 This program controls printing of the maneuver out- 
put. When a printout point is reached., subsidiary variables such 
as q^ and q are calculated.. 

7. INFOL This is a routine to perform linear interpolation 
for a single variable. 

8 . LAGRA This program carries out quadratic interpolation for 
nonequispaced points. The function value and/or its first deri- 
vative is produced. 

9. AREA This program controls the calculation of the rays, 

area, and. pressure. It saves the ground intersection points and. 
provides for the calculation of all 0 's . 

10. DERIZ The calculation of the integrands for the ray and. 
area integrals is carried, out by this routine. 

11. OUTPU This program calculates the age function using 
Simpson’s rule for integration. It provides for printing ray, 
area, and age data at printout points. 
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Figure 13. Program organization 
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12 • DERPR This program calculates the integrand required by 
the pressure calculation. 

13. OUTPR Evaluation of nhase distortions £ 

-functions, 

and pressure data, and printing of them, are carried out by this 
routine . 

14. NODE This routine provides for the solution of a set of 
n first-order ordinary differential equations or integrals. A 
Runge-Kutta-Gill (RKG) starting procedure and a highly stable 
Adams-type predictor-corrector method are used. 

15. PASS This is an auxiliary routine used by NODE when 
evaluating integrals. 

16. GROUN This program retrieves the ground intersection data 
and reshuffles it to provide tables of ray location and maximum 
pressure in. order of increasing time for each maneuver point. It 
then determines the ground intersection contours for equal time 
on the ground. 

17. INTEP This program provides for the simultaneous linear 
interpolation of x , y , and pressure for a given value of time. 

18. MOVE This routine moves tabular data from location J to 
location K . 


19. RESTA This is a routine used to save data for restarting 
the problem at a later time. 


2 °. INITAL This routine provides for starting a problem from 
previously saved data. 


Three sense switches (SS) may be used, as follows: 
Sense Switch Function 


SS-1 Restart solution from data saved on Tape 9 . 

SS-5 Print out during all integrations if trouble 

is suspected. 


SS -6 Stop calculations at end of next ray trace 

and pressure calculation. Save data on Tape 9 
for future restart. 


The names of variables (FORTRAN symbols) which are used, in 
the program are identified in Table 1 . Comments which will be 
listed by the program when appropriate are shown in Table 2, in- 
cluding suggestions for correcting errors. 
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Program Utilization and Instructions 

This section describes how to set-up the inputs for utili- 
sation of the digital program and explains applications of the 
equations presented in THEORETICAL ANALYSIS. Details which are 
needed for setting up a card, input declare presented in Tahl s 
3 4 and 5 and explained in the following paragraphs. .Table 3 

shows how input data are arranged on punched cards, giving the 
data formats and notes relating to their set up. Table pre- 
sents definitions of the input nomenclature and option numbers, 
whereas Table 5 lists the units used in the three choices of in- 
put units. The following subsections are labeled to correspond 
to related sections of THEORETICAL ANALYSIS. 

The Atmosphere .- The inputs of tabular atmospheric data are 
made on the sets of cards labeled Cards 5* 6, and 7 in 3. 

The temperature and either the pressure or the density are tabu- 
lated for various arbitrary altitudes in ascending order. Either 
the pressure or the density table is used, not both A separate 
card is used for each altitude. The units are seiected according 
to Table 5, choosing a number for IIJNIT according to Table 4. If 
the 1962 U.S. standard atmosphere is desired, ISTAN - 1, and 
Cards 5 and 6 are omitted from the deck. The wind magnitude and 
wind direction are listed for various altitudes on Cards J. .The 
altitude on Cards 5, 6, and 7 may be geopotential or geometric, 
selecting a cor responding digit for IALT (Table 4) . 

lated T Snear e iy P wSreal ' tET^nf tabSe fs^n^LpSate/quaSatJ^ly. 

able errors in the interpolation. 

Aircraft Maneuvers. - Physical data regarding the ^aircraft 
are placed on Card“ The altitude HG is. the height of the 
o-round above sea level whereas HO is the initial altitude of th 
aircraft above sea level. The initial values of the x,y coordi- 
nates of the aircraft are set to zero within the program. Its 
position and altitude at later times are calculated Th^wing^ 
tions (2) based on the maneuver inputs of Cards 8. The wing 
loading WS is used only in an auxiliary calculation of force co- 
Pfficilnts Cr and Ct - C D (eqs. 3 and 4), and does not enter 
into the sonic boom calculation here (except indirectly through 
-hhe F-function values). The aircraft length LA is. used as. a 
reference length, whereas the radius of the earth is used in con- 
verSoSs^ betSfen' geopotential and geometric altitudes m the at- 
mospheric functions. 

The maneuver information is placed on Cards 8 as functions 
of increasing flight time. For Option 1 (IM = 1) the axial and 
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lifting load, factors, Mach number, flight path angle , heading 
angle , and bank angle are input using a separate card for each 
flight time. The inputs should, of course, be made consistent j 
that is , the Mach number and angles should be values which corres- 
pond to the load factors n T and n L . Equations (5) and (7) 
are used to obtain the time derivatives. For Option 2 (IM = 2 ) 
the load, factors are not required. These tables are interpolated 
quadratically, so that at least three maneuver cards are needed. 
The ^ time derivatives are obtained by calculating slopes of quad- 
ratic curves fitting the data. On Card 9 are given the initial 
time and. subsequent time intervals along the flight path for 
which rays are to be initiated. 

The force coefficients are determined using equations (3) and 
( 4 ). with the dynamic pressure calculated in the program. In 

Option 2, the load factors used in these equations are given by 
n L ~ cos 1 and n T - sin 7 to provide approximations to com- 
plete the digital program output listing. This substitution is 
exact only for nonaccelerating, straight flight. 

Initial Wave Normals . - The ray parameters are specified by 
the time ta (Card 9) and the azimuth angle 0 . This angle is 
first set to zero within the program, and. then incremented, by an 
amount A 0 (DELFI) after each ray-path integration until 0 MfiY 
(FIMA.X) is exceeded, as input on Card. 10 . The negative incre- 
ments —A 0 are then taken until ^MAX d® exceeded.. For a prob- 
lem which is symmetric in 0 , the input DELFI can be negative 

to avoid, duplicate calculations of the pressure signatures at 

+0 and -0 . The values of v and sin 0 O are evaluated 

using equations (9) and (10) . 

Mach Conoids and Ground Intersections . - Each time a ray 
trajectory intersects the ground, its location, arrival time, and 
peak pressure are saved in a special table. This peak pressure 
is selected from the array of pressures (Ap as functions of 
phase i) which are determined, but not saved, in other parts of 
the program.!} After all maneuver, ray tracing and pressure 
calculations have been completed, these saved data are interpo- 
lated to obtain intersection data at specific arrival times at 
the ground, as follows. 

1 ) These pressure values may be useful indicators of the 
actual pressure increment at each location, but they are not 
accurate. The peak values have not yet been adjusted., at this 
point of the procedures, for the distortion of the signature and 
the development of shock-wave structure . The user must judge the 
value of these interpolated, peak pressure listings and, for 
accurate values, properly process the signature data given in 
the detailed tables for each ray-tube calculation. 
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The program selects from the table the minimum arrival time 
t for the set of rays which left the aircraft at the second 
maneuver time (t a +i) • Then the program proceeds to the previous 
ray set corresponding to t a and interpolates linearly to find 
the ray locations and peak pressures for the same arrival time t . 
There are generally two of these, corresponding to starboard . and 
port ray azimuths. The program then finds a new minimum arrival 
time for the set of rays corresponding to t a+2 and interpolates 
the preceding sets (t a and t a+1 ) . This process continues 

through each set of rays corresponding to subsequent aircraft 
maneuver times. 


An example of output for this ground intersection table is 
shown in figure l4. These data were obtained for large ray azi- 
muth increments of 20° and, therefore, yield somewhat crude inter- 
polations. The maneuver flight times at which ray calculations 
were initiated were at one-second increments. The solid-line 
curves show ground intersection locations for sets of rays for 
these maneuver times. The dashed-line curves show ground inter- 
section interpolations for the fixed arrival times of 37.6, 38.4, 
39.3, and 40.1 seconds. 


Snell's Law and Ray Tracing. - The ray tracing initial inte- 
gral ioiTTtep^zMEP7 _ lirdr^rin : rTntervals (ZPRN) for the solution 
’ * > i8) are specified on Card 10. Equation 16 is cal- 

* M I • I _ ^ 1 4 "1 n TA /A rl 


of equations 
culated for 


A quantity c 


os 


= c 


are listed at the head of 


+ u n o is also obtained.. 


_ - . " Q • ^ ~OQ 

The values of 0 , v , c OQ and c 0 

each ray tracing table (Appendix B) . The wind components u n 
and ut are obtained using equations (15) and cos 0 is ob- 
tained using equation (17) • The ray trajectory integration. step 
is halved or doubled as the calculations proceed from the air- 
craft to the ground, according to predictor-corrector error 
limits and comparisons made in sub-routine NODE. .Listings of ray, 
area and age variables are made at descending altitude increments 
specified by ZPRN. These listings are also given at the altitude 
corresponding to the next integration interval below the selected 
altitude called PRIN input on Card 10 „ 


The calculations along a given ray are stopped at the ground 
(z = 0), or if the ray-tube area converges to zero, or if the ray 
angle decreases to within 2° of horizontal (cos 0 — 0.9994) . The 
additional displacement during the final 2° to the true horizon- 
tal point is approximately 


57.3°/rad d(a - u n )/d(-z) 


with the quantities evaluated at the horizontal point. In cer- 
tain special problems, more than one Ar contribution may be 
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Figure 14. Typical ground intersection data 
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needed. For standard atmosphere at sea level without winds, the 
correction is 


A? = 1.93 stat mi (3.1 km) 


Rav-tube Area. - The partial derivatives of__c Q and v 
W ith and 0 are calculated in subroutine AREA using equa- 
tions (27). (28) and. (29). Variables dependent on altitude z 

and the integrals (eqs. (24)) are evaluated using subroutines 
DERIZ and NODE, and the various terms of the area equation (do) 
are combined in subroutine OUTPU. The values of the ray-tube 
area are listed in the output table as a function of altitude z . 

Flow N ear the Aircraft . - The program input of the aircraft 
F -function is made of Cards 4. Fq is tabulated as a function 
of length parameter L/LA . The input format permits specifying 
the Fi vs L/LA function for several azimuth angles 0 r (Phi) 
relative to the aircraft vertical (butt-line). The program se- 
lects that Fq vs. L/LA curve given in the input table which has 
0 r closest to the current computed value 0 r = 9 - 9a j where 

0 is the ray parameter and 0 a is the bank angle of . the 
craft. The factor Ff is input on Card 10, and multiplies Fq 
before it is listed in the output tables (eqs. (34)). 

This format for inputing the aircraft’s signal requires . that 
the dependency of F on 0 r , C^ , weight and Mach number is 
determined independently of this sonic boom program. Other sub- 
routines for calculating variations in Fq .may be substituted., 
however, by a programmer. Also, the dimension specifications m 
this program may need to be changed to allow for larger input 
tables . 

This program calculates pressure signature parameters at 
each given value of L/LA and at three additional intermediate 
stations, using linear interpolation to determine the intermediate 
values of the F-function. These extra data are. useful for de- 
termining the pressure curve and the ?/- curve which yields the 
shock locations. 

Geometr ic Acoustics and Blokhintsev Invariance .— The quantity 
y r 45 ) ) is evaluated in subroutine DERP R . The phase 4 

(eq. (4i)) is determined in OUTPR. These are listed in the output 
tables (Appendix B) for corresponding values of L/LA . 


Signal retor tion and Age Variable .- .Equation (46) for. the 
age variable must be rewritten for the. digital computations in 
order to properly handle its initial singularity. . / P 

singularity occurs because the initial ray-tube area is zero. Its 

evaluation is as follows: 
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The age variable is 



d(-z) 


(64) 


where z Q is the airplane altitude where the ray begins and z-, 

is any lower altitude. G through G-jo represent combinations 

o ^ ^ he area ^^ion (26) which are used in the program 

and defined in Table 1. & 


Let 


Q(z) = 


G5G4 


/ 


z - z 
o 


a/A 


(65) 


so that 


T 


- 1.2 


a c 


1 Q(z q ) dz 1 

— + 

V~z - z 


Q(z) - Q(z q ) 


z o ' ~° “ ~ Z D ^ Z ° " Z 


dz 


iy 2 [-2Q(z 0 )v^-r^ + J 


Q(z) - Q(z o ) 


z o 


dz 


(66) 


This correctly yields t & = 0 at z = z Q when Q(z a ) is deter- 
mined In the vicinity of z , let A = A-, (z n - z). 2 ) Then from 
equation ( 26 ) and Table 1 , 10 


A l (z o) “ a o G l G 3 G 8 G 12 V + G 0 G 2 G 6 Gi 


7 


(67) 


Then Q(z q ) can be evaluated as 


«■«>> “ 7^ 


(68) 


S g 5 g 4 ^ 


2) 


The area A is linear in z for a uniform atmosphere. 
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in the digital program, the age is calculated ih subroutine OUTPU 
from equation (66) using equations (65), (67), an “ t 1 ' , fhA lth 

9imn^on integration procedure for the integral. X . 

the P ray-tube area and ray coordinates as a function of altitude z . 

The phase variable ^ is calculated at the ground (z = 0) 
in subroutine OUTPR and is listed m the output table. 

Shock Location.- The integral S Q (|) (eq.(54)) ^ he iraliioq 

function o/Uq) (55)) are obtained in subroutine OUTPR. Vela 

of these variables are listed in the output table at corresponding 


>1 


where SINT is 


K- 


input on 


o(^) 

Card 10, 


values of £ and 
The reflection factor 
the pressure, also in subroutine OUTPR. ine 
increment Ap (eq. (57)) and pressure ratio 


is jsf(li) 
to multiply 
pressure 
(where p 
respectively, in 


and S 
is used 
resulting 
Pi = Ap/p 


is ambient pressure) are listed as DELP and PI, 
this output table. 

As described in THEORETICAL ANALYSIS, Shock Location the user 
must locate the shocks and the shock jumps by plotting XUi) 

Ap versus given in this output table . 


Conversion to Other Computers 


The program listing of Appendix B is the program used for 
computations on the CDC-6600. The following remarks are pertinent 
to conversion to other computers: 

1 Input and output units are standard FORTRAN units 5 and 
6. For other units, change KUNIT = 6 and LUNIT = 6 in 
subroutine INPU. If a different working tape than tape 
unit 8 is used (saves ground intersection data), change 
MUNIT = 8 in INPU. 

2 Some computers use FORMAT statements requiring the holli- 
rith notation 1 , whereas the symbol * is used here for 
the CDC-6600. 


o To fit within the core storage limitations of . the IBM- 
1130, links, locals and socals have been required . 

4. Monitor control cards must be added to be compatible 
with the compiler. 


5. 


6 . 


The sense switch subroutine SSWTCH should be checked for 

acceptability . 


'he dimension statements pertaining to the tabular input 
ata may need to be changed to accommodate larger tables. 
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COMPUTATION RESULTS 
Sample Printout 

Appendix B presents a reproduction of a complete printout 
for a typical solution. Various sections of the printout a?e 
delineated on the right-hand side by letters A through E? 

, . ^ Se n' tion .^. is bhe in Put data which is listed for identifica- 
tion and verification. The number of significant figures in the 
is ting might be less than those on the input data cards and 
stored within the computer. These inputs ?epreJent an aSplane 
flying m an easterly direction O = 90 °') at 20 non r - h ^ 

ft ?^ lne z t ”V ea b Maeh U 5 - ^VoSd'altXS 

angles^ 1 '- rf 1 F "ri U !i2 10n tabulated ion two ray azimuth 

thf aircraft The p 4 ? 5 f ° r 19 stab ions along the length of 

tine aircraft. The F-function table differs for the two 

oniy at station L/LA = 0 . 60 . The wind direction A from the west 

Special uAits t ttnwTt nd) i N bUt fV zero s P eecU The inputs are in 
Stffdftef f fi, f f ” V’ atmosphere is the 1962 D.S. 

biandard (ISTAN - 1). Rays are to be initiated at the aircraft 

ne-second intervals of flight time (DELTA T PRINT) Data are 

20° te terff f ° r - ayS ln , ltlated at azimuth angle intervals of 
. , to a maximum value of ± 6 0° (MAXIMUM PHI) using an 

tDElTA 1 7F te f rat i°S lx } t f rval of approximately 300 ft (91 m) 

> uTA Z). Ray-tube information is requested at 2000 ft (610 m'l 
intervals of altitude (PRINT INTERVAL) and at the next inte^ra- 
tron interval below 3500 ft ( 1.1 km) (PRINT OUT POINT) The 

SffsT a S wte° n °L K P ” hiCh “Itiplies tee p^lure tecre- 
(eq. (^Ij.) 1 ) whereas the F-function multiplier Fp is 3.784 

,, Section B illustrates Maneuver Data output. These include 
the state variables of the aircraft at the time t Q when the rav 
tracing is started. Force coefficients, wind components and 
ynamic pressure are _ also listed. Maneuver Option 1 was 5 selected 

makin°- uh of tv, ? ri Y a ^ 1V !f S ^ 3 ^ 3 ^ were calculated 

makirit, use of the load factor data in the input section (n T n L ) . 

Section C represents the data listed during propagation of 

fro ?. th ? aircraft to the ground. Thf firffhf idL- 
tifies the particular azimuth angle 0 and. the values of v, 

c oo and c Q at the aircraft altitude. Below this is a table of 
ffai tliul°e n ’ z . ime> ray ‘ tUbe area, age, and cos 9 as fnjtifnf 

Section D is the pressure signal data at the ground. Data 
are given for each integration interval of the V’-curve. It is 
important to have a large set of points so that the pressure and 
curves can be plotted accurately. In this table, ? F is the 
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V (L) is the signal invariant 


input value F* multiplied by F± , yt'T' is the phase t n , 
r^n (Ur')'), XI o is the phase \ , XI 1 is rne pause ti , 

£ e ?‘ VJ integral of V F U) , and s is the value of €he 
columns . 


Both 

variable 


Ap and are to be plotted against the independent 

if as shown in figure 15 for a M = 1.2 solution. 

=> 1 




The shock locations are defined wherever J ° ros ^ s V^uppermost 
proceeding along the curve with “creasing on the upp 

branches °yt° u ndary (sup ^^xa^le office if) because,/ 
abscissa (i = -0.004 m tne & nressure 

otherwise appear to be multivalued. 

K-SS'^SSS ?alu^ e ofX d l^Sr^d n peak 

pressure for constant arrival time at the ground. 
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Typical Results 


T 7 Un*^u VarlOUS solutlons have been calculated to check the equations 
which were programmed and to compare a few sample results with 
previous analyses. Two characteristics which San bl evaluated 
analytically to provide checks on many factors and terms in the 
equations, but excluding effects of winds, are 

a) wKSh^ra SreJd C i? CUla h Path (ref - 38) ln an atmosphere 

icn has a speed of sound varying linearly in z 

b) ff ea Wh l Ch 13 com P ute h for an atmosphere with a 
onstant temperature will increase quadratically in the 

lateral distance parameter (ax 2 + bx + c if the air 

SS a ”e^v t ?^S n ft SOUt !? ) - More ° ver > lf the airerert Is 

in steady flight (time derivatives all zero), the area 
will increase proportionately to the decrease in z 

representative°examples? rately ^ theSe anal ^ lc ^ lts ^r 

been SlS?ated n to S ^ tl0 ? S f or typical practical situations have 
Appendix r 1 21?° d ! monstrate the program using the F-function of 

li P ?Sr Steady ?ught at ^ fig ?f eS 15 3113 

directly bSlS; thS airereft® U 1 th ? gr ° Und lntersectl ° n 
for m = l p hnth 9 = °1* In com Panson with the data 

the greater aee b of%hP m " C ^ rv ? an f the Pressure curve reflect 
ulo f rea ? er of the M = 3 signature, as the pressure profile 

on bright fOT^ I" ha" aVe Meed, the two shocks 

e risnc lor M = 1.2 have merged into only one at m = 3 

, , Results using an atmosphere which has a larger lapse rate 

at Vooo fWq’Tk f tem P er y ure of 9°°P at the "ground and -75°F 

at 80 000 ft (9.1 km), are shown in figure 17 for several Mach 
numbers. This atmosphere is designated ATM A -3 in reference 48 
e P re ^sure ratio here is the maximum pressure listed for the 

TJSZ? :lT P ° h Te e . re diVlded by the mLlra ™ P«-BU re i f n r a the 


Tailwind effect 
18 for some initial 
figure 19 (see ref. 
pressure listed for 
ponding value in a s 
give pressure ratios 
a more serious sonic 


s on signal strength are summarized in figure 
calculations. The wind is the mean zonal of 
f? e . Pressure ^ ratio here is the maximum 
c ailwind solution divided by the corres— 
tandard no-wind atmosphere. A headwind would 
larger than 1.0 near Mach 1.2 and represents 
boom environment. 
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a) Pressure signature 



b ) ^-curve 


Figure 16. 


Distorted, pressure signature and 
standard atmosphere, M = 3 




-curve at ground 
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100 



0 0.02 0.04 

Wind speed, km/sec 

Figure 19. 

Mean zonal wind 


The preceding results were for rays 
starting directly beneath the aircraft, 
0=0. For increasing 0 , the maximum 
lateral extent of ray intersections with 
the ground is shown in figure 20 for a 
standard atmosphere. The effect of wind 
on lateral range is shown in figure 21. It 
is seen that lateral range is considerably 
reduced by a headwind. Results with a 
nonstandard atmosphere are shown in figure 
22 for atmosphere B-5 of reference 48 
(multiple temperature inversions in tropo- 
sphere ) . 

The wind and temperature effects 
brought out by these examples on sonic 
boom intensity and maximum lateral range 
are similar to those presented in more 
detail in reference 48. Specifically, 
only small differences from standard atmo- 
sphere, no wind solutions are indicated at 
aircraft speeds larger than M = 1.3 • 

The results of the sample problem 
listed in Appendix B provide an example of 
a solution with an airplane maneuver. For 
an axial acceleration of 0.2 g's and a 
lift acceleration of 0.8 g's in a turn, 
the maximum pressure is 6l^ greater than 
for nonaccelerating straight flight at 
Mach 1.5. Thus, while atmospheric varia- 
tions may be expected to be essentially 
unimportant to sonic boom intensity at 
speeds above M = 1.3 ^ maneuver effects 
can continue to have a strong effect. 


CONCLUDING REMARKS 

This report has presented an analysis of sonic boom propaga- 
tions in a stratified atmosphere with winds, and a computer program 
based on the analysis. The analysis and the program take into 
account maneuvers of the aircraft and yield actual pressure signa- 
tures without the common N-wave approximations. Some additional 
theoretical discussion and historical notes have been included to 
make the report more useful as an exposition of the fundamental 
theory. Sample computer input and output data are presented with 
results of some preliminary calculations made using the program. 
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Figure 22. Lateral range with atmosphere B-5 


The digital program was written to allow a number of input 
options and to be readily adaptable to various computers. Complete 
definitions and instructions for its use are presented. 

Preliminary results have been shown which indicate the capa- 
bilities of this new program. It can be applied to analyzing . the 
effects on sonic boom strength of changes in atmospheric conditions, 
wind profiles, and aircraft maneuvers. It is, therefore, useful to 
the aircraft designer to show effects of aircraft geometry and 
initial signature, to the environmental science services (meteorol- 
ogists) to show significance to atmospheric variations for specific 
aircraft, and to the governmental aviation authorities and the 
airlines operations analysts to provide a basis for specifying 
permissible flight maneuvers and flight profiles within sonic boom 
overpressure constraints. 

It is possible and desirable, of course, to make certain 
additions and alterations in the future to the program. One of 
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these may be to provide equations at the beginning of the program 
for directly calculating the aircraft F-function, now required as 
a known input. Another is to extend the automatic calculations to * 
provide directly the /-curve intersections (shock locations) and 
shock strengths , and to provide automatic plotting of the final 
pressure signatures. These are now manual operations. The present" 
program, in common with other known techniques, cannot provide 
pressure information when a caustic develops, nor can it account 
for atmospheric turbulence. These are two important areas 
requiring further research. 
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TABLE 1. 


COMMON 

A 

AO 

ACC "1 
ACCUR > 
PHSIG J 

AGE 

ALT 

ARA 

BETA 

BYHAL 

REDUC 

C 

CN 

CONI 

C0N2 

CON3 

CON4 

COO 

CO 

CT 

CTH 

CU 

D 


NAMES OP VARIABLES IN PROGRAM 


Speed of sound 

A at the initial point of ray 

Accuracy criteria for integrals evaluated using NODE 
Age variable t 

Altitudes in altitude/temperature input table 
Area of ray tube 

Halving suppression switch for integration 
Coefficients for Simpson's evaluation of AGE 

COS T) 

y G5O~cos 2 0 o /tan p 1 also use a as coefficients in 

, n | AGE calculation 

sm p/c Q -cos 6 j 

AGE/2 

(G50/G4 • ARA) 1//2 

a /cos 0^ 
o / o 

c - u 
oo no 

cos 0 
cos 6 

o 

cos JJL 

Atmospheric density for a given altitude z 



DADZ 


DEL 

DELPI 

DELP 

DELPT 

DELZ1 

DPMAX 

DSIM 

DTPRN 

EETA 

ENDNO 

ENDVA 1 
ENDVL J 

ETA 

ETADZ 

P 

FA 

PIA 

FIMAX 

PI 

FLAG 

GAMR 

GDOT 

GR 

GO 


da/dz 

Current step in 0 
Angle increment for 0 

Maximum and initial step size for ray-area calculation 

Print interval for ray-area calculation 

Maximum initial step size for ray-area calculation 

Maximum value of pressure for all F's at a specified 0 

Step for Simpson's rule integration for AGE 

Print interval for maneuver calculation 

Wind direction for a given z 

Number of points in integration range 

Value of independent variable at end of integral 

Wind direction in input table 
dr]/dz 

Two-dimensional array of F-function input 
Bank angle at initial point of ray 
Bank angle in input table 
Maximum value of 0 

Ray azimuth angle 0 measured from vertical plane 
Set to -1 to terminate integration 
Flight path angle 7 in input table 
dy/ dt a 

Value of 7 for a given time 
cos 9/sin 0 
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G1 

G2 

G3 

G 4 

G5 

G50 

G6 

G7 

g8 

G9 

G10 

Gil 

G12 

G13 

H 

HNODE 

HG 

HW 

HO 

HI, H2, 
HI 


G 0 3 /a 2 


1 + (sin 7/sin [i sin 0 ) . 

. ^ , ° ) also used as temporaries 

cos 7 cos | a sin 0/ cos 0 Q J 


sin 9 cos 0 
pa 2 


p a 
H o o 


2 


o 


cos [i/cos 0 
cos p sin 7 + sin p cos 7 cos 0 
sin 0 /cos-* 0 

2 

cos p sin 0 Q sin 0/cos 0 Q 


sin 0 Q sin 0 


2 


cos 7 sin 0/cos 0 Q 

cos p cos 7 sin 0 

sin p cos 7 + cos p sin 7 cos 0 

Step size used for integration 

Height of ground above sea level 
Height for input wind table 

Initial value of aircraft height above sea level 

H3 H4 H5, H7: Variables representing terms and factors 

' in equation for ray-tube area: 

d.0 _ 


= u 


to dt 


'a 


a o G 8 G 7 u to G lly dt a 


a G Q G no - u, G. 


o 


8 13 " to 9 / dt 


dX. 


'a 


- (7 q sin 7 r / sin p^ - a^/cos 6 q 
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H2 


QPL. . n T d U + r r d 7 
dt a + G ll dt a + G 6 G 10 dt^ 


H3 “ " u to G 6 G 7 " a o G 8 G i2 


H4 


= - g 6 g 7 


H5 


- G 3 V + u tQ G 2 


H7 

I 

IALT 

IM 

I OUT 

IPRIN 

IRH 

IRHO 

I STAN 

I TEMP 

IUNIT 

KTIME 

L 

LAMDO 

LAMR 

LA 


(H1)(H4) - (H2)(H3) 

Counter for Simpson's rule and various DO loops 
Altitude geopotential/geometric switch 
Maneuver option 1 or 2 
Switch for first z printout 
Switch for printing at a selected z 

Switch to indicate if ray is horizontal or area is zero 

Altitude , temperature , pressure/density switch 

Standard/tabular atmosphere switch 

Type of temperature units input 

English/metric/special units switch 

Number of ground intersections stored 

DO instruction parameter for steps in S-integral 

Value of ip for a given value of time 

Input table of heading angles (f) 

Length of airplane 
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LUNIT 


M 

ML 

MM 

MTIME 

MUDOT 

MUNIT 

MU 

N 

NNODE 

NALTS 

NEW 

NIM 

NL 

NODUB 

NODUM 

NPHI 

NSAVE 

NT 

NTAU 

NWIND 

P 

PHI 

PRESS 


Output unit number 

Mach number input table 

Position in PHI table for F-function 

Mach number for a given time 

Switch to indicate if first maneuver step 

dp/dt a 

I/O unit on which ground data are saved 
Mach angle p. 

Number of integrals to be evaluated simultaneously 

Number of entries in the altitude/temp/density table 

Heading angle v of wave propagation vector 

Number of entries in maneuver table 

Load factor input table, lift direction 

Switch to suppress doubling step during integration 

Switch to suppress printing of bad points during 
integration 

Number of 0 1 s for the F— function table 

Number of auxiliary variables calculated during 
integration 

Load factor input table, axial direction 
Number of F-functions for each 0 
Number of entries in wind table 
Pressure for a given z 

Input table of 0 r 1 s for which F-values are tabulated 
Pressure Input table 
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PRIN 


PSIDT 

Re 

RE 

RHO 

RHOO 

SN 

ST 

STH 

SU 

SUM 

T ], 

TNODE J 

TAU 

TEMP 

TIME 

TINE) 

TINIT 

TLAS 

TPRIN 

TP 

TT 

UNO 

UTO 


SGlsctsd value of z at which additional printing 
is desired 

df/ dt & 

Radius of the earth 

Reflection factor 

Atmospheric density input table 

Density at the initial point of ray 

sin v 

sin 6 

sin 0 

o 

sin p 

Used in collection of terms for AGE calculation 

Used in integration step size determination for NODE 

Length ratio L/LA in input table 
Temperature input table 
Time at initial point of ray 

Time at current integration point in maneuver 
Time at which maneuver is to be started 
Last value of time at which interpolation was made 
Next value of time for a maneuver printout 
Temporary location 

Input table of times for maneuver table 

Component of wind speed in direction of wave propagation 
Component of wind speed transverse to wave propagation 
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V 

VAR 

WW 

VWX 

vwxx 

VWY 

VWYY 

V¥ 

WS 

XM 

YM 

Z 

ZINIT 

ZPRIN 

ZSIMP 


Speed of the airplane 

Two-dimensional array used to save back points in 
integration 

Input table of wind speeds 
X-component of wind at maneuver point 
X-component of wind 

Y-component of wind at maneuver point 
Y-component of wind 
Value of wind speed at a given z 
Wing loading 

x at initial point of ray (a maneuver point) 

y at initial point of ray 

z at current point of ray 

z at initial point of ray 

Next value of z to be printed 

Value of z for next Simpson's step 


INPU 

FF Factor to rescale F table 

I DO instruction parameter 

J DO instruction parameter 

KUNIT I/O unit for input 

ATMO 

D Density at height z 
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H z plus height of ground 

HGLTB Array for base heights 

HGP H in table value units (geopotential meters) 

I DO Instruction parameter 

J Layer number being used 

K1 Array of coefficients 

K23 Array of coefficients 

P Pressure at height z 

PB Array of P-g coefficients 

RHOB Array of p-g coefficients 

TB Array of tg, coefficients 

TEM1 Difference between HGP and table value 

TEM2 Quantity used in calculation of D 

TH Temperature at height z 

THDZ dTH/dz 

VS Speed of sound at height z 

Z Current height above ground 

ANGLE 

ARRAY Array to be checked 

I DO instruction parameter 

MAX Number of entries in ARRAY table 

MAXS MAX - 1 
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DERM2 


BLANK 

COSGA 

COSLA 

J 

SI NLA 

TH 

TM 

0UTM2 

BLANK 

CL 

COSDI 

COSGA 

COSMU 

Cl 

NLZ 

NTZ 

Q 

SINDI 

SINGA 

TERM1 

TH 

THDZ 

TM 


Dummy variable for LAGRA call 

cos 7 

cos 

Sense switch value 
sin i/ 

Temperature from standard atmosphere program 
Temperature at a given z 


Dummy variable for LAGRA call 

Lift coefficient 

cos (f - T)) 

cos 7 

cos |T 

Axial force coefficient 
Value of NL at point z 
Value of NT at point z 
Dynamic pressure 
sin (f - T) ) 
sin 7 

Subterm in derivative calculations 

Temperature at z from standard atmosphere program 
dTH/dz by quadratic interpolation 
Temperature for a given z value 


103 



VDOT dV/dt a 

VG Velocity of aircraft relative to ground 

INPOL 

FT Interpolation value 

I DO instruction parameter 

K Place in table where found 

T Point at which interpolation value is needed 

X Table containing point at which interpolation occurs 

Y Table to be interpolated in 


LAGRA 

ANSI 


ANS2 


A1 

A2 

A3 


} 


B1 

B2 

B3 


! 


Cl 

C2 

C3 


} 


L 


N 

NN 

NSWIT 


Interpolated value 
Interpolated derivative value 

Differences between input T and tabular values 
in table TIME 

Differences between three consecutive values of T 
in table TIME 

Terms used in evaluation formulas 

DO instruction parameter 
Number of entries in each table 
N - 1 

-1 only derivative, 0 both, 1 only function 
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T 

TABLE 

TIME 

AREA 
COSET 
CTH2 
I DEL 
SINET 

DERIZ 

AST 

J 

W 

OUTPU 

TAUA 

PER PR 
J 

OUTPR 

PI 

S 


Point at which interpolated value is to be found 
Table to be interpolated in 
Table which contains values of T 


COS T] 

cos 0 C 


o 


Number of print steps 
sin rj 


1 


a • sin 6 

Setting of sense switch 5 

Value of wind velocity at given z 


Integrand for AGE integral 


Setting of sense switch 5 


Ap/p 

SINT - 0.5 x age v|{ |) 
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SINT 


XIO 

XII 


NODE 



CK 


COMPD 

COMPE 

COMPT 

COMPY 

I,J,K 

NN 

NSWHF 

PERR 

R 

TEMP 

TEMPA 

XSAVE 

GROUN 

PI 

I 


fV E U) 

Linear phase (time) | 
Actual phase (time) 


Coefficients for RKG starting procedure 
Subterm in RKG evaluation 

Name of subroutine which calculates integrands 

Name of subroutine accessed before doubling and halving 

Name of subroutine to produce output 

Used in differential equations 

DO range parameters 

Number of integrals plus saved quantities 

Rehalving switch 

Maximum error at step 

Subterm in RKG evaluation 

Predictor value of variable 
th 

Error for n equation 

Value for independent variable when step was halved 


0 for given ground data 

position in tables at which current data are stored 
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IC 

IMIN 
I SAVE 

J 

JMAX 

K 

KERR 

KI 

KK 

KMAX 

KN 

LOOPS 

N 

NMNSM 

NPHSW 

NREC 

NRMAX 

P 

PVAL 

T 

TMIN 

TVAL 

X 

XVAL 


Step for moving and for looking up in tables 

Position in main table where current table starts 

Place in tables where values are stored during 
permutation 

Place in current T-table which has minimum value of time 

Place in tables where current table ends 

Place current value is stored in N-table 

Indicates when there is no TVAL in table 

Next position in table and DO parameter 

Table number 

Last table used for current curve 
Starting place of current table 
Place moved from 
Table of table lengths 

Switch to indicate first negative 0 for maneuver point 
Switch to indicate first maneuver point 
DO range parameter 

Number of ground intersections saved 

Table of pressure values at the ground 

t InL 

Pressure value at time TVAL in KK table 

Table of time values at the ground 

Minimum value of time for the current table 

Time along curve being plotted 

Table of x-values at the ground 

t h 

x-value at time TVAL in KK table 
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Y 

YVAL 

INTEP 

IN 

INF 

IS 

KK 

N 

NERR 

P 

PANS 

PERCT 

T 

TVAL 

X 

XANS 

Y 

YANS 

MOVE 

J 

K 


Table of y-values at the ground 

f-v, 

y-value at time TVAL in KK table 


Place in T-table where current table starts 
Place in T-table where current table ends 
DO instruction parameter 

Place in table N which gives position in table T 
Table of table lengths 

Set to 1 if no time TVAL in table , otherwise 0 
Table of pressure values at the ground 
Pressure at time TVAL 

Distance between table entries for values 
Table of time values at the ground 

Entry in time table at which interpolation is made 
Table of x-values at the ground 
x at time TVAL 

Table of y-values at the ground 
y at time TVAL 


Position in tables moved from 
Position in tables moved to 
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RESTA 


C 

INITAL 

C 


COMMON block 


COMMON block 


QRU¥ih;%L IS 

OF POOP Ql.’A?..:rs 
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TABLE 2. 


COMMENTS GENERATED IN PROGRAM 

Interpretation 

Altitude in tables must increase 
monotonically; correct data and 
rerun the problem. 

Time in maneuver table must 
increase monotonically; correct 
data and rerun problem. 

Maximum 0 must be zero if 60 
is zero; correct data and rerun 
problem. 

Tables should have at least three 
entries. Revise data and rerun 
problem. 

Linear interpolation is attempted 
outside the range of the table. 
Execution continues but answers 
will not be correct. Revise 
table or the problem statements. 

6. AREA ZERO Area in ray tube diminishes to 

zero. Current ray calculation 
is terminated. Age may not be 
exactly correct. 

7. RAY HORIZONTAL Ray has turned to become within 

2° of horizontal. Current ray 
calculation is terminated. 

8. NEGATIVE SQUARE ROOT IN Current ray calculation is 

AGE terminated. 

9. ALTITUDE LARGER THAN Altitude outside range of standard 

52 000 METERS atmosphere program; execution is 

terminated . 


Comment 

1. ALTITUDE VALUES NOT 
INCREASING 


2. TIME IN TABLE NOT 
INCREASING 


3. DELTA PHI IS ZERO WITH 
PHIMAX NOT ZERO 


4. TABLES TOO SHORT 


5. ABSCISSA VALUE OUTSIDE 
RANGE 
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TABLE 3. 


FORMATS FOR INPUT DATA 
Card 1 Title card 


Card 2 

(All input 

is fixed point. 

right 

justified) 


Cols 1-5 

6-10 

11-15 

16-20 

21-25 

26-30 

31-35 

IUNIT 

I STAN 

IALT 

I TEMP 

IRHO 

IM 

NPHI 

Cols 36-40 

41-45 

46-50 

51-55 




NTAU 

NALTS 

NWIND 

NIM 




Card 3 

(All input 
is E20 . 7 } 

is floating point F10 
right justified) 

. 3 , except 

last 

Cols 1-10 

11-20 

21-30 

31-40 

41-60 



ws 

LA 

HG 

HO 

RE 





Cols 1-10 11-20 21-30 

PHI L/LA F 

List L/LA and F for a given 0 r , then proceed to next larger 0 r 
and list corresponding L/LA and F . 

If loading nonstandard atmosphere,, use Cards 5 or 6; otherwise 
omit . 


Cards 5 to (4 + NALTS ) ( F10 . 3 ) 

Cols 1-10 11-20 21-30 

ALT TEMP PRESS 



Cols 1-10 11-20 21-30 

ALT TEMP RHO 
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Cards 7 to (6 + NWIND) ( F10 . 3) 

Cols 1-10 11-20 21-50 

HW WW ETA 

Cards 8 to ( 7 + NIM) (F10.5) 

If IM = 1 


Cols 1-10 

11-20 

21-30 

31-40 

4i-50 

51-60 

61-70 

T 

NT 

NL 

M 

GR 

PSIR 
= LAMR 

FIA 

If IM = 2 







Cols 1-10 

11-20 

21-30 

31- 

-40 

51-60 


T 

M 

GR 

PSIR = 

= LAMP 

FIA 


Card 9 (F10. 

3) 






Cols 1-10 

11-20 






TINIT 

TSTEP 






Card 10 (F10.3) 






Cols 1-10 

11-20 

21-30 

3i~4o 

41-50 

51-60 

61-70 

ZSTEP 

ZPRN 

DELFI 

FIMAX 

PRIN 

RF 

FF 


NOTES : 

1. Altitude means geometric altitude except in input when 
geopotential option is used. 

2. Angles in input and output are in degrees. 

3. Output listings are English units. The wind speed in output 
listings has units ft/sec, whereas in input the units are knots. 
The unit for pressure is lb/ft^ . 

4. Pressure, density, temperature and F-function tables are inter- 
polated linearly, whereas other input tables are interpreted 
quadratically . Therefore the tables require a minimum of 
three input cards. The inputs must be selected to avoid wrong 
quadratic representation between the specified data. 

5. The F-function used in the computation is the tabular data for 
the 0 nearest the 0 being computed. 
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6. For a steady-state flight, a minimum of three maneuver cards 

must be used with arbitrary time T increasing from the first 
card to the last. To avoid duplicate solutions, select TSTEP 
so that TINIT + TSTEP is larger than the latest time T . 

7 if both A 0 and 0 max are input as zero, only the 0=0 ray 
Is computed. If A 0 is input as zero and 0 max 1S n ™ zer °’ 
an error statement is given and no solution is made. Ii A 0 
is input as a negative number, only the rays for 0=0 and 
negative 0's are computed. 


113 



TABLE 4. 

INPUT DEFINITIONS 


IUNIT : 

ISTAN: 

IALT: 

I TEMP: 


IKHO : 

IM : 

NPHI : 
NTAU : 
NALTS : 
NWTND : 
NIM : 
WS : 
LA : 
HG : 
HO : 
RE : 

PHI : 


-1- -ENGLISH 

0- -METRIC Input units; outputs are English units 

+1- -SPECIAL 

1 - -1962 US Standard. Atmosphere 

2- -Tahular atmosphere 

1 — Geopotential altitude 1 input for tabular atmosphere 

2 — Geometric altitude J and wind data 

1- -Fahrenheit 

2- -Centigrade 

3- -Rankine 

1- -Altitude, temperature, pressure table 

2- -Altitude, tempera ture, density table 

1- -Maneuver, option 1 

2 - -Maneuver option 2 

Number of PHI 1 s ( 0 ) in F-function table 

Number of length parameters L/LA in F-function table 

Number of entries in IRHO table (Cards 5 or 6 ) 

Number of entries in wind table (Cards 7 ) 

Number of entries in maneuver table (Cards 8 ) 

Wing loading, W/S 
Length of aircraft 

Altitude of ground above sea level 

Altitude of aircraft above sea level 

Radius of earth (e.g., 2.089007 x 10^ ft or 
6.367293 x 106 m) 

Ray azimuth angle relative to aircraft z-axis, 0 r 
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FI : 
L/LA : 

F : 

ALT : 
TEMP : 
PRESS : 
RHO : 

HW : 
VWW : 
ETA : 

T : 

NT : 
NL : 
M : 
GR : 
PSIR : 

FIA : 
TINIT : 

TSTEP : 

ZSTEP: 
ZPRN : 

DELFI : 


Ray azimuth angle relative to vertical plane, 0 
Stations at which F-function is specified 
Aircraft signature function 

Atmospheric temperature and either pressure or 
density as functions of altitude above sea level. 
Altitude (ALT) may be geopotential or geometric. 

Wind speed and direction from which wind is coming as 
functions of altitude (HW) above sea level. Direction 
is the heading angle measured eastward from north. 

Aircraft flight time 

Axial load factor 

Lift (normal) load factor 

Mach number 

Flight path angle above horizontal'] relative to the 

_ . r atmosphere. 

Heading angle measured clockwise j includ i ng w j. n d 
from north (eastward) 

Bank angle of aircraft, 0 a 

Initial time at which first ray tracing is to be 
calculated (must correspond to any T input) 

Time interval along aircraft flight path at which 
ray calculations will be initiated 

Initial ray-area integration interval of time (DELZl) 

Altitude intervals at which ray data are to be listed 
(DELPT) 

Increments of FI (ray azimuth angles) at which 
calculations are desired 


FIMAX 

PRIN 

RF 

FF 


Maximum FI desired 

A selected altitude at which ray data are desired 
Ground reflection factor (normally 1.8 to 2.0) 
F-function factor to multiply input parameter F ± 
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TABLE 5. 

INPUT UNITS AND CONVERSION FACTORS 



OPTION (- 1 ) 
English 

OPTION (0) 
Metric 

OPTION (+ 1 ) 
Special 

Atmospheric temper- 
ature, T 

. °F or °C 
: or °R 

O-n, o„ 

F or C 

or °R 

o n o„ 

F or C 

or °R 

Atmospheric pressure 

, p: lb/ft 2 

N/m 2 

mbars 

Atmospheric density. 

p :slug/ft^ 

kg/iP 

O 

kg/irn 

Wind speed, V or u 
w 

: knots 

m/ sec 

knots 

Wind direction, r\ 

: deg 

deg 

deg 

Length or altitude 

: ft 

m 

ft 

Wing loading, W/S 

: lb/ft 2 

kg/m 2 

lb/ft 2 




CONVERSION FACTORS 

ft to m 


: multiply by 0.3048 

lb/ft 2 to 

mbars 

: 0.47880258 


N/m 

r\ 

: 47.880258 


j 

kg/m 

: 4.88242 

slug s/ft 3 

to kg/m^ 

: 515.379 

knots to 

ft/ sec 

: I. 6878 I 


m/sec 

: 0.514444 

0 

ct 

O 

O 

o 


: °C = (5/9) (°F - 32) 

°F to °R 


: °R = °F + 459.67 
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APPENDIX A 


PROGRAM LISTING 


PROGRAM SONIC 

EXTERNAL DERM2 .OUTM2 .PASS 

REAL M.LAMDO .MM ( 2 5 ) . LAMR t 2 5 ) . NT ( 25 ) . NL { 25 ) . LA . MU . MUDQT 
DIMENSION PH I ( 20 ) *TAU <20)»F(20.20).FIA<25). 

1 ALT < 25 ) *T EMP < 25 > . PRESS (25)»HW(25).VVW(25)*ETA(25)*TT<25>.GR (25) 
COMMON NNODE .HNODE .TNODE.ACC .NODUM.BYHAL » NODUB * ENDNO * ENDVL * NSAVE * 
5FLAG.TIND.VAR< 14.8) 

COMMON I UN IT . I STAN . I ALT . I TEMP . I RHO. IM.NPHI.NTAU.NALTS.NWIND.NIM 
1 .WS.HG.HO.RE.PHI .TAU.F.FIA.ALT .TEMP. PRESS. RF * HW . VVW . ETA .MU .MUDOT » 
2TT»MM,GR»LAMR.NT*NL*LA»TINIT .GDOT .DTPRN.DELZl .DELPT . F I MAX . DELF I . PR 
3 IN .MT I ME .KT I ME »M . V . VWX . VWY . GAMR .PS I DT » RHOO » LAMDO . TPR I N . T IME.ZINIT* 
5XM.YM.A.BETA. LUN IT .DADZ.VW.EETA »VWDZ * ETADZ * MUNIT .FI .FA.DEL.TLAS 

EXECUTIVE PROGRAM FOR THE SONIC BOOM CALCULATION 
CALL SSWTCHI 1 .NNODE ) 

GO TO (4.1 ) .NNODE 
C THE REAL BEGINNING OF A DATA CASE 

1 CALL I NPU 

2 F I =0 » 

DEL = DELFI 

C MAIN PROGRAM FOR MANEUVERS 

80 I F ( T PR I N-TT ( N I M ) ) 70.70.20 
70 IF < ABS< TPRIN-TT (NIM) >-.01 ) 20.20.10 

20 NNODE = KTIME+1 

WRITE(MUNIT) F I . F I . F I . F I . F I 

NODUM = MUNIT 

NODUB = LUN IT 

REWIND MUNIT 

CALL GROUN 

GO TO 1 

10 TIND = TPR I N 
TLAS = -9999. 

VAR (1.1) = XM 
VAR (1.2) = YM 
VAR (1*3) = ZINIT 
GO TO (30.40) .MT I ME 
30 MT I ME = 2 
ENDNO = 0. 

CALL DERM2 
CALL 0UTM2 

C CALL RAY TRACING THE FIRST TIME 


GO TO 

3 

40 TPR I N 

= TPR I N + DTPRN 

NSAVE 

= 0 

ENDVL 

= TPR I N 

NNODE 

= 3 

FLAG = 

0* 

ACC = 

1. 

TNODE 

= X. 

NODUM 

= 0 

NODUB 

= 0 

BYHAL 

* *5 

ENDNO 

= 10* 


50 CALL NODE ( DERM2 .DERM2 .0UTM2 .PASS) 

3 CALL AREA 
GO TO 2 

4 CALL INITAL 
GO TO 3 
END 
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SUBROUTINE INPU 

REAL M.lAmM .MM (25) .LAMR125) . NT ( 25 ) . NL ( 2 5 ) .LA.MU.MUDOT 
DIMENSION PHI ( 20 ) *TAU ( 20) , F( 20 *20 ) .FI A ( 25 ) * 

1 ALT (25) * TEMP ( 25 > »PRESS(25) . RHO( 25 ) .HW125 ) * 

2 VVW ( 2 5 ) ♦ ET A ( 2 5 ) , T T ( 25 ) » GR ( 2 5 ) 

COMMON N.H.T ..ACCUR , NODUM , REDUC. NODUB , ENDNO. ENDVL * NSAVE ♦ FLAG* Z* VAR 
1(14*8 ) 

COMMON I UNIT . I ST AN . I ALT ♦ I TEMP ♦ I RHO » IM.NPHI . NT AU *NALTS ♦ NW I ND * N I M 
1 *WS*HG,HO*RE»PHI » T AU . F . F I A » AL T . TEMP ♦ PRESS ♦ RF * HW . VVW . ETA . MU .MUDOT » 
2TT»MM*GR*LAMR*NT *NLtLA*TINIT *GD0T *DTPRN *DELZ1 .DELPT *F I MAX *DELF I .PR 

3 IN »MT I ME *KT I ME * M * V *VWX * VWY *GAMR . PS I DT *RHOO » LAMDO.TPR IN.TIME.ZINIT* 
5XM* YM.A .BETA.LUNIT .DADZ » VW » EETA . VWDZ . ETADZ »MUN I T 

EQUIVALENCE (PRESS! 1 ) *RHO( 1 ) ) 

C 

C 

C I UN I T = NEGATIVE NUMBER FOR INPUT VALUES IN ENGLISH SYSTEM 
C LUNIT = 0 FOR INPUT VALUES IN METRIC SYSTEM 

C IUNIT = POSITIVE NUMBER FOR INPUT VALUES IN ESSA 
C 

C ISTAN = 1 FOR STANDARD ATMOSPHERIC DATA 
C ISTAN = 2 FOR TABULAR ATMOSPHERIC INPUT 
C 

C I ALT = 1 FOR ALTITUDE IN GEOPOTENTIAL UNITS 
C I ALT = 2 FOR ALTITUDE IN GEOMETRIC UNITS 
C 

C ITEMP = 1 FOR TEMPERATURE INPUT IN DEGREES FARENHE I GHT 
C ITEMP = 2 FOR TEMPERATURE IN DEGREES CENTIGRADE 
C I TEMP= 3 FOR TEMPERATURE IN DEGREES RANKINE 

C 

C I RHO =1 FOR ALT I TUDE .TEMPERATURE* PRESSURE TABLE 
C I RHO = 2 FOR ALTITUDE, TEMPERATURE. DENSITY. TABLE 
C 

C IM = 1 FOR MANEUVER OPTION 1 , TABLES OF LOADFACTORS. GAM.PSI AND BANK 

C IM IM = 2 FOR MANEUVER TWO. FOR MACH. GAMMA AND PSI TABLES 

C 

C NIM IS THE NUMBER OF ENTRIES IN THE MANUEVER INPUT TABLES 
C 

C NPH I IS THE NUMBER OF PHI'S IN THE F FUNCTION TABLE 

C NTAU IS THE NUMBER OF TAU ' S IN THE F FUNCTION TABLE 

C 

C NALTS IS THE NUMBER OF ENTRIES IN THE ALT I TUDE . TEMPERATURE ♦ 

C PRESSURE/DENSITY TABLE 
C 

C NWIND IS THE NUMBER OF ENTRIES IN THE WIND VELOCITY TABLE 
C 

C RE IS THE RADIUS OF THE EARTH 

C 

C 

C*****READ CONTROL CARD 
KUNIT = 5 
LUNIT = 6 
MUNIT = 8 
READ (KUNIT. 4) 

READ ( KUN I T . 1 ) 

1 IUNIT. ISTAN. I ALT .ITEMP. I RHO. IM.NPHI .NTAU . NALTS ,NW I ND » 

2NIM 

1 FORMAT ( 11 (2X. 13) ) 

C*****PRINT CONTROL VALUES 
WR I TE ( LUNI T » 8 ) 

WRITE ( LUN I T .4 ) 

WRITE (LUNIT. 2) 
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WR ITE ( LUNIT ) 

1 ( ALT ( I ) *TFMP ( I ) * RHO ( I ) * I =1 ,NALTS ) 

702 WR ITE ( LUNIT * 15 ) 

WRITE ( LUNIT t 16 ) 

1 ( HW ( I ) * VVW ( I) *ETA( I ) *1=1 *NW IND ) 

GO TO ' 70,60 ) * I V; 

60 WR ITE ( LUNIT. t61 ) 

WRITE (LUNIT *62 ) 

1 ( T T ( I ) * MM ( I ) , GR ( I ) , L AMR ( I ) , F I A ( I ) * I = 1 , N I M ) 

GO TO 900 

70 WRITE (LUNIT* 71) 

WR I TE ( LUNIT * 72 ) 

1 ( TT ( I ) * NT ( I ) * NL ( I ) * vm< I ) * GR ( I ) *LA^R( I ) * F I A ( I ) * I = 1 * N I R, ) 

900 WR ITE ( LUNIT *91 ) T I N I T » DTPRN 

WR I TE ( LUN I T * 92 ) DEL7 1 * DELPT * PR I N * OELF I * F I WAX *RF *FF 
IF(DELPI) 160* 155*160 

155 IF(FIMAX) 158*156,158 

158 WRITE (LUNIT *159) 

159 FORMAT ( 2 0 X * 39HDELT A PHI IS ZERO WITH PHI WAX NOT ZERO) 

CALL EXIT 

750 WPITF (LUNIT *251) 

251 FORMAT ( //20X30HATLITUDE VALUES NOT INCREASING ) 

CALL FXIT 

757 WRITE (LUNIT *253) 

753 FORMAT ( //70X *28HT IME IN TABLE NOT INCREASING ) 

CALL EXIT 

156 DELFI =-.00005 
C 

C 

C CONVERSIONS 

r 

C 

G CONVERSION TO GEOMETRIC UNITS IF ALTITUDE IN GEOPOTENTIAL UNITS 
C 

160 GO TO (700*204) * ISTAN 

204 I P ( N A L T S— 2 ) 150,150*205 

205 GO TO( 201*2 0 3 ) * I ALT 
201 DO 181 I = IfNALTS 

181 ALT(I) = ( ALT( I )*RE) / (RE-ALT( I ) ) 

CONVERSION TO 7 

203 DO 183 1=1 *N ALTS 
183 ALT ( I ) = A LT ( I ) -HG 

CHECK FOR TABLES TO HAVE INCREASING VALUES 
DO 206 I = 2 * NAL T S 

I F ( AL T ( I ) - ALT ( I - 1 ) ) 250*250,706 

206 CONTINUE 

200 GO TO ( 1820 * 1821 )* IALT 

1820 DO 182 I=1*NWIND 

182 HW ( I ) = ( H W ( I ) *RE ) / ( RE-HW ( I ) ) 

1821 DO 184 1 = 1 * N W I N D 

1R4 HW ( I ) = HW(I) -HG 

DO 20 7 I = 2 * N W I N D 

IF(HW(I)- H W ( I — 1 ) ) 250 *250,207 

207 CONTINUE 

DO 208 I = 2 * N I M 

IF (TT ( I )-TT ( 1-1 ) ) 252 *252 *208 

208 CONTINUE 

CONVERSION TO FNGLISH SYSTEM IF NECESSARY 
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1 IUNIT.ISTAN.IALT.ITEMP.IRHO, IM.NPHI .NTAOtNALTS.NWIND.NIM 

IF1NWIMD-2) 150*150.151 

150 WRITE! LUNIT » 152 ) 

CALL EXIT 

152 FORMAT ( 20X » 20HTABLES TOO SHORT ) 

151 I F ( N I m -2) 150.150.154 

C READ WING LOADING. LENGTH OF AIRPLANE. ALTITUDE 

C OF THE GROUND. ALTITUDE OF THE AIRCRAFT. AND RADIUS OF THE EARTH 
154 READ ( KUNI T . 3 ) 

1 WS.LA.HG.HO.RE 

3 FORMAT(4F10.3.E20.7) 

8 FORMAT ( 1H1 ) 

4 FORMAT (80H 

1 > 

'C READ AIRCRAFT SIGNATURE. F FUNCTION 
RE AD ( KUN I T . 5 ) 

1 ( ( PHI ( J ) .TAU! I).F(I.J).I=1 »NTAU ) . J=1 .NPHI ) 

GO TO! 103.22 ) . I STAN 

C READ ALT ITUDE. TEMPERATURE ♦ PRESSURE /DENS I TY TABLE 
22 GO TO (101.102) . I RHO 

101 READ ( KUNI T . 5 ) 

1 (ALT! I ) .TEMP! I ) .PRESS! I ) . 1 = 1. NALTS) 

GO TO 103 

102 READ ( KUN I T . 5 ) 

1 (ALT! I ) .TEMP! I) .RHO! I ) . 1=1 .NALTS) 

C REAR WIND VELOCITY TABLE 

C HW IS THE ALTITUDE. WW IS THE WIND SPEED. AND ETA IS DEGREES FROM 
C NORTH 

103 READ ( KUNIT ♦ 5 ) 

1 (HW! I) . V V W ( I).ETA(I).I=1.NWIND) 

CALL ANGLE(ETA.NWIND) 

IF(IM-l) 7.7,6 

6 READ ( KUNIT. 55 ) 

1 ( TT ( I ) .MM! I ) ,GR ( I ) » L AMR (I).FIA(I)»I=1»NIM) 

GO TO 90 

7 READ ( KUN IT » 57 ) 

1 (TT! I ) .NT! I ) »NL( I ) »MM( I ) ,GR( I ) . LAMR ! I).FIA(I).I = 1.NIM) 

CALL ANGLE ( F I A .N IM ) 

90 READ ( KUN I T » 5 ) TINIT.DTPRN 
CALL ANGLE (GR.NIM) 

CALL ANGLE ( LAMR.NIM) 

READ(KUNIT»58)DELZ1.DELPT .DELFI .F I MAX . PR I N . RF , FF 
58 FORMAT ( 7F10. 3 ) 

FI MAX = APS ( F I MAX ) 


♦****dr int OUT ALL 


INPUT DATA 


WR ITE ! LUNIT » 141 ) WS.LA 
WRITEILUNIT.14) HG.HO.RE 
WR ITE! LIJN I T . 1 3 ) 

WRITE (LUNIT. 12) 

1 ( (PHI ( J) . T AU ( I ) » F ( I » J ) » 1 = 1 .NTAU) . J=1 .NPHI ) 

GO TO ( 202.222 ) , ISTAN 
222 GO TO! 111,112) , I RHO 
111 WR ITE1LUNIT .10) 

WRITE! LUNI T .9 ) 

1 (ALT(I), T EM P(I), PRESS !I), 1=1, NALTS) 


GO TO 202 

112 WRITE (LUNIT .11 ) 
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c 


I F < IUNIT13030.301.338 


CONVERSION FOR METRIC TO ENGLISH 
301 HG =HG/ . 3048 
H0=H0 / • 3048 
LA=LA/ • 3048 
WS = WS/4. 88242 
RE=RE / < 3048 
DO 381 1 = 1* NW I ND 

HW ( I ) = HWt I ) /.3048 

381 VVW ( I ) = ( VVW ( I ) / • 3048 ) 

GO T 0 ( 405 * 406 ) . I RHO 

405 DO 483 I = l.NALTS 
ALT ( I >=ALT< I 1/.3048 

483 PRESS! I ) “PRESS t I )/ 47 *88026 
GO TO 333 

406 DO 484 1=1 *NALTS 
ALT! I ) “ALT ( I 1/.3048 

484 RHO! I 1 =RHOl I 1 /515.379 
GO TO 333 

CONVERSION FROM ESSA 

338 GO TO (3030*3338) »ISTAN 
3338 GO T0(601*602)*l RHO 

601 DO 681 1 = 1 * NALTS 

681 PRESS! I (“PRESS! I )/ .4788026 
GO TO 3030 

602 DO 682 I=1*NALTS 

682 RHO! I )=RHO( I 1/515.379 

C CONVERT WIND SPEED FROM KNOTS TO FEET PER SECOND 
3030 DO 782 I=1*NWIND 

782 VVW! I) “VVW! I 1*1.688944 

C CONVERT TEMPERATURE TO RANKINE 

333 GO TO ( 303*3031 1 ♦ ISTAN 
3031 GO TO (901.902.303) *ITEMP 
901 DO 981 1=1. NALTS 
981 TEMP! I 1 “TEMP ( 11+459.7 
GO TO 303 

902 DO 982 I=1*NALTS 

982 TEMP! I ) = <9./5.)*<TEMP( I 1+273.151 
C CONVERT ALL DEGREES TO RADIANS 
303 DO 781 1=1 . NPH I 
781 PH 1(1) “PHI < I )*.1745329E-01 
DO 783 1=1 . NW I ND 

783 ETA! I) “ ETA ( I 1 *. 1745329E-01 

784 DO 786 I=1*NIM 

GR( I 1 = GR ( I 1 *• 1745329E— 01 

786 LAMR! I)=LAMR< I 1 *. 1745329E-01 

787 DO 788 I =1.NIM 

788 F I A ( I 1 =FIA( I 1 * . 1 745329E-0 1 

2 FORMAT ( /10X.6HI UN I T“ * 13* 10X * 6H I STAN= * I 3 * 1 OX * 6H I ALT“ I3.10X. 

1 6HITEMP=.I3*10X* 5HIRHO=* 13* 10X*3HIM=. 13.// 

2 10X * 6HNPH I = * I 3 * 10X *6HNTAU= » 1 3 » 10X » 6HNALTS“ » I 3 * 10X * 6HNW I ND= * I 3 * 

3 10X. 5HN I M= .13 ///) 

5 FORMAT (3F10.3) 

9 FORMAT! 17X.F7.0.22X.F5.0.15X.E20.7) 

10 FORMAT!// . 16X » 8HALT ITUDE.19X.1 1HTEMPERATURE ♦ 2 IX ♦ 8HPRESSURE 1 

11 FORMAT!// . 16X .8HALTI TUDE. 19X . 1 1HTEMPERATURE « 22 X » 7HDENS I TY 1 

12 FORMAT! 20X.F5.2 » 1 OX . F5 . 2 . 10X . F 1 0. 5 1 
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13 F0RMAT(21X* 3HPHI *12X»4HL/LA >16X.1HF) 

14 FORMAT ( 10X»*HG=*>F7.0*10X»*HO=*.F7.0»10X»*RE=*E20.7/// ) 

15 FORMATt// 19X*8HALTITUDE*8X*10HWIND SPEED ♦ 5X * 9HD I RECT I ON ) 

16 FORMAT (20X*F7.0*8X*F8.2*9X * F6 • 1 ) 

55 FORMAT ( 5F 1 0 • 3 ) 

57 FORMAT ( 7F 10 • 3 ) 

61 FORMAT ( / / /23X*4HTIME*17X»7HMACH NO * 14X * 5HGAMMA * 14X * 4H PSI*10X*10HB 
1 ANK ANGLE) 

62 FORMAT (20X»F10.3»10X»F10.3»10X.F10.3*10X»F10.3»10X»F10.3) 

71 FORMAT ( // 13X.4HTIME*15X.*NT*.15X*NL*10X**MACH NO*12X*GAMMA*l IX 

1 »*PSI*10X,*BANK, ANGLE*) 

72 FORMAT ( 7F17. 3 ) 

91 FORMAT { //1X*13H1NIT1AL T I ME = . F 10 . 3 * 1 3X . 14HDELT A T PR I NT* . F14. 5 / / ) 

92 FORMAT ( 1X*8HDELTA Z= . F 1 4. 5 . 1 4X ♦ 15HPR I NT I NTERVAL* . F 10 . 3 » // * 1 X 

2 16HPRINT OUT PO I NT* » F 1 0 . 3 * 10X * 5HDPH I = * F 10 . 3 * 1 OX * 1 2HMAX IMUM PH I = * 

3 F10.3//1X.19HREFLECTION FACTOR =, F10.3. 7X.10HF FACTOR =,F10.3) 
141 FORMAT ( 10X . 13HWING LOAD I NG= * F 10 . 3 » 1 OX * 16HA I RPLANE LENGTH* ♦ F10 . 3 

1///) 

785 TPRIN = TINIT 

DO 789 J = 1 * NPH I 
DO 789 I = l.NTAU 
789 F ( I » J ) = FF*F ( I ♦ J ) 

ZINIT = HO-HG 
XM = 0. 

YM = 0. 

MTIME * 1 
KT I ME = 0 
REWIND MUNIT 
RETURN 
END 
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SUBROUTINE A TMO ( Z »T H »P » D ♦ VS ) 

r!!!l M*LAMD0^»MM(25) » L AMR ( 2 5) » NT ( 2 5 ) * NL ( 2 5 ) .LA.MU.MUDOT 
DIMENSION HGLTB( 6 ) tP* ( 5 > *RHOB ( 5 ) tTP ( 5 ) 

DIMENSION PH I ( 20 ) *TAU ( 20 ) t F ( 20 » 20 ) * F I A ( 2 5 ) *RHO ( 25 t , 

1 ALT(25)*TEMP(2^)*PRESS(25)fHW(25)*VVW(25)*ETA(25) tTT(25)tGR(25) 

1 CO'MMON NNODE .HNODE.TNODE.ACC » NODUM ♦ BYHAL » NODUB * ENONO • ENDVL * NS AVfc ♦ 

l.WS.HG.HO.RE.PHI.TAU.F.FIA.ALT.TEKP.PRESS.RF* HW . VVW . ET A .MU . MUDOT . 
2 IT *MM , GR * LA MR *NT *NL .LA *T I N I T » GDOT .DTPRN.DELZ1 » DELPT » F I MAX iDELF I ♦ R 

3 I N *MT I ME * KT I ME >M * V »VWX « V'WY » GAMR *PS I DT » RHOO * LAMDO »TPR I N * T I ME * Z I N I T * 

5XM , YM * A t RET A * LUN I T »DADZ .VW.EETA » VWDZ * ETADZ *MUN I T 

EQUIVALENCE! PRESS( 1 ) *RHO( 1 ) ) 


C H IS THE HEIGHT IN METERS FOR WHICH A CORRESPONDING TEMPERATURE. 

C PRESSURE AMD DENSITY WILL RE COMPUTED 

C HGLTR»K1 tK23*TR*PBtRHOR ARE STANDARD ATMOSPHERE DATA ARRAYS 

C T H » P * D * V S * ARE THE TEMPERATURE . PRESSURE t DENS I TY * AND SPEED OF SOUND 

C TO RE COMPUTED AND RETURNED TO THE CALLING PROGRAM 
C 
C 

c 


HGLTB! 1 ) 

= C • 

K 1 < 1 ) 

= -. 225571E-04 

K 2 3 ( 1 ) 

=-5.255894 

TB ( 1 ) 

=518.670 

PB ( 1 ) 

=2116.21695 

RHOB ( 1 ) 

=2 . 37692E— 03 

HGLTB! 2 ) 

=11000. 

K 1 ( 2 ) 

= 0. 

K.2 3 ( 2 ) 

=. 157689E-03 

TP ( 2 ) 

=389.970 

PR ( 2 ) 

=472.68 

RHOB! 2 ) 

=7 . 06 126E— 4 

HGLTR ( 3 ) 

=20000. 

<1 ( 3 ) 

= . 46 1 574E— 0 5 

K23 ( 3 ) 

=34.163426 

TR ( 3 ) 

=389.970 

PB t 3 ) 

=114.346 

RHOB ( 3 ) 

=. 170817E-3 

HGLTB ( 4 ) 

=32000. 

K1 (4 ) 

= . 12245E-04 

K2 3 ( 4 ) * 

12.2012? 

TB ( 4 ) 

=411.57 

PR (4 ) 

=18.129 

RHOB l 4 ) 

= 2 . 56609E— 0 5 

HGLTB! 5 ) 

=47000. 

K 1 ( 5 ) 

= 0. 

K 2 3 ( 5 ) = 

. 12622733E— 03 

TB ( 5 ) 

=487.17 

PB ( 5 ) 

= 2.3163 

RHOR ( 5 ) 

=2 . 76983E— 6 

HGLTB (6) 

=52000. 

FIND LAYER NUMBER 


H = Z + H G 

HGP=0.304R*H/(1.0+0.3048*H/6356766.0) 

IF (HGP-53000. > 1.5.5 
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1 DO 20 1*1,5 

1 F ( HGLTB ( I )-HGP>2»3«7 

2 J= I 

20 CONTINUE 
GO TO 7 

3 J=I 

GO TO 7 

5 WR I TE ( LUN I T ♦ 6 ) 

6 FORMAT ( //20X33HALTITUDE LARGER THAN 52000 METERS) 

CALL EXIT 

C 

C TO COMPUTE TEMPERATURE, DENSITY, PRESSURE AND SPEED OF SOUND 
C 

7 TEM1 =HGP-HGLTB ( J ) 

TH=TB ( J ) 

THDZ = TH*K 1 ( J ) 

I F ( K 1 ( J ) >103,104,103 
103 TEM1=1.0+K1 ( J)*TEM1 
TH=TH*TEM1 

TEM2 = TEM1** ( -1 .0-K23 ( J ) ) 

D= RHOB ( J ) *TEM2 
P =PB( J)*TEM1*TEM2 
GO TO 105 

00104 TEM1 =EXP ( -K23 ( J ) *TEM1 ) 

P = PB ( J ) *TEM1 

D =RHOB( J)*TEM1 
105 VS=49.02057*TH**0.5 

DADZ = ( 1201.5081*THDZ)/VS 

RETURN 

END 


SUBROUTINE ANGLE (ARRAY >MAX ) 

DIMENSION ARRAY (10) 

MAXS = MAX - 1 
DO 4 I = 1 , MAXS 

9 IF(ABS(ARRAY( 1+1) -ARRAY! I ) 1-180.) 4,4,8 
8 IF(ARRAY(I+1) -ARRAY ( I ) ) 10,4,11 

10 ARRAY ( I +1 ) = ARRAY ( I +1 ) +360. 

GO TO 9 

11 ARRAY ( I + 1 ) = ARRAY { I + 1 ) -360. 

GO TO 9 

4 CONTINUE 
RETURN 
END 
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SUBROUTINE DERM2 

REAL M. LAMDO .MM(25) .LAMR! 25) . NT t 25 ) * NL ( 2 5 ) *LA*MU*MUDOT 
DIMENSION PH I ( 20 ) * T AU ( 20 ) ,F<20.20>.FIA(25> .RHO< 25) ♦ 

1 ALT ( 25 ) *TEMP (25)*PRESS<25)»HW(25).VVW(25)»ETA(25)»TT<25)»GR(25> 
COMMON NNODE »HNODE . TNODE . ACC . NODUM . BYHAL * NODUB .ENDNO.ENDVL *NSAVE» 
5FLAG.TIND.VAR! 14.8) 

COMMON IUNIT .ISTAN.IALT .ITEMP.I RHO. IM.NPHI .NTAU.NALTS.NWIND.NIM 
l.WS.HG.HO.RE.PHI »T AU .F.FIA.ALT * TEMP .PRESS. RF . HW . V VW , ET A , MU .MUDOT . 
2TT.MM.GR. LAMR. NT. NL.LA.TINIT.GD0T.DTPRN.DELZ1.DELPT.FIMAX.DELFI.PR 
3IN.MTIME.KTIME.M.V .VWX .VWY . GAMR .PS I DT . RHOO . L AMDO . TPR IN.TIME.ZINIT. 
5XM.YM.A.BETA.LUNIT .DADZ.VW.EETA.VWDZ.ETADZ. MUNIT .FI .FA.DEL.TLAS 
EQUIVALENCE! PRESS! 1 ) ,RHO( 1 ) ) 

3 IF (TIND-TLAS ) 5.60.5 

5 CALL LAGRAITIND.NIM. TT.GR. GAMR. BLANK. 1 ) 

CALL LAGRAtTIND.NIM.TT.LAMR.LAMDO. BLANK. 1 ) 

CALL LAGRA ( T I ND . N I M . T T . MM » M » BLANK . 1 ) 

TLAS = TIND 
SINLA = SIN! LAMDO ) 

COSLA = COS (LAMDO) 

COSGA = COS (GAMR) 

60 CALL LAGRAtVAR! 1.3) .NWIND.HW.VVW.VW . BLANK. 1) 

CALL LAGRA! VAR! 1 .3 ) . NW I ND . HW , E T A , EE T A . BLANK . 1 ) 

GO TO (40.50). ISTAN 
40 CALL ATMO ( VAR ( 1 . 3 ) . TH . P » D . A ) 

GO TO 55 

50 CALL I NPOL ( ALT .TEMP .VAR(1»3) .TM) 

A = 49.020576*SQRT (TM) 

55 V = A*M 

VWX = VW*S I N ( EETA ) 

VWY = VW *COS ( EETA ) 

VAR (8.1) = V *COSGA*SI NLA-VWX 
VAR (8. 2) = V *COSGA*COSLA-VWY 

VAR (8.3) = V *S I N ( GAMR ) 

CALL SSWTCH ( 5 » J ) 

GO TO ( 100.200) »J 

100 WRITE ( LUN IT .155) T I ND . VAR ( 1 . 1 ) . VAR ( 1 . 2 ) . VAR ( 1 . 3 ) 

155 FORMAT!* TIME = * E15.8 .5X.*X =* » E 1 5 . 8 ♦ 5 X . *Y =* . E 1 5 . 8 » 5X » *Z =*E15 
1.8) 

200 RETURN 
END 
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SUBROUTINE QUTM2 

REAL M*LAMDO ♦ MM ( 2 5 ) * L AMR ( 2 5 ) * N T ( 2 3 ) ♦ NL ( 2 5 ) * L A 
1 *MSG*MUDQT*MU*NTZ*NLZ * MDOT 
DIMENSION PHI (20) * T AU ( 2 0 ) * F ( 2 0 * 20 ) * F I A ( 2 5 ) *RHO ( 25 ) • 

1 ALT ( 25 ) *TEMP ( 25 ) *PRESS ( 25 ) *HW ( 25 ) *VVW (25)*ETA(25)*TT(25)*GR(25) 
COMMON NNODE *HNODE ,T NODE * ACC* NODUM*BYHAL * NODUB * ENDNO * ENDVL * NSAVE • 
5FLAG*TIND*VAR ( 14*8 ) 

COMMON I UN I T * I STAN* I ALT-* I TEMP . I RHO* IM*NPH I *NT AU * NALTS * NW I ND *N I M 

1 * WS *HG *H0 *RE * PHI * TAU * F * FI A * ALT * TEMP *PRE5S *RF * HW *VVW * ETA *MU*MUDOT* 

2 T T *MM*GR * LAMR * NT * NL * L A * T I N I T * GDOT *DTPRN*DELZ1 *DEL.PT *F I MAX .DELF I * PR 

3 IN *MT I ME *KT I ME *M * V *VWX * VWY *GAMR *PSI DT * RHOO * LAMDO * TPR IN*TIME*ZINIT* 
5 X M * Y M * A * B E T A * L UN I T * DADZ * VW * EE T A * VWDZ » ET ADZ ♦ MUN I T * F I *FA*DEL*TLAS 

EQUIVALENCE ( PRESS ( 1 ) *RHO< 1 ) ) 

IF (ENDNO) 20*20*100 
?0 FLAG = -1. 

F!H = VAR (1*3) +HG 
GO TO ( 11*13) * IM 

11 WRITE(LUNIT*7) 

G = R E / ( RE + HH ) 

G = 3 2.257 2 4*G*G 

GO TO 12 

13 WRITE (LUNIT * 1 ) 

10 CALL LAGRA(TIND*NIM*TT*GR*BLANK*GD0T*-1 ) 

CALL LAGRA (TIND*NIM*TT * LAMR *BLANK*PSIDT*-1 ) 

CALL LAGRA ( T I ND * N I M ♦ T T * MM * M ♦ MDOT * -1 ) 

12 WRITE (LUNIT *2) T I ND * HH 
M S Q = M*M 

BETA = SORT ( MSQ-1 . ) 

MU = AT AN ( 1 • /BETA ) 

G 1 = 57.29578*GAMR 

G2 = 57.29578*LAMD0 
G3 = 5 7.2957 8*MU 

WR I TE ( LUN IT*3)VAR(1*1J *VAR(1*2)*VAR(1*3) »G1*G2*G3 
GO TO (50*55) * ISTAN 
50 CALL A T M 0 ( VAR ( 1 * 3 ) * TH * P *RHOO * A ) 

GO TO 60 

55 GO TO (70*71)* IRHO 

71 CALL INPOL ( ALT *RHO * V AR ( 1 * 3 ) * RHOO ) 

GO TO 72 

70 CALL INPOL ( ALT*PRESS»VAR( 1 *3) *P ) 

CALL INPOL(ALT*TEMP*VAR(l*3)*TM) 

RHOO = P / ( 1 7 1 6 • *T M ) 

7 2 CALL LAGRA ( VAR ( 1 * 3 ) * NAL TS * ALT * TEMP * BLANK * THDZ *-l ) 

DADZ = ( 1201 .5081*THDZ ) /A 

60 Q = .5*RHOO*V*V 
COSMU = COS (MU) 

SINGA = SIN(GAMR) 

COSGA = COS(GAMR) 

CALL LAGRA ( VAR ( 1 * 3 ) *NW I ND * HW * VVW* BLANK * VWDZ *~1 ) 

CALL LAGRA ( VAR ( 1 * 3 ) * N W I ND * HW * E TA * B LANK * E T ADZ *-l ) 

CALL LAGRA ( T I ND * N I M * T T * F I A * FA * BLANK *1 ) 

GO TO. (61*62) *IM 

62 VDOT = A*MDOT+(MSQ*DADZ*SINGA)*A 
Cl = (WS*SINGA)/Q 
CL = ( WS*COSGA)/Q 

GO TO 63 

6 1 CALL LAGRA (TIND*NIM»TT » NL » NLZ * B LANK * 1 ) 

CALL LAGRA ( TIND*NIM*TT*NT*NTZ*BLANK*1) 

COSDI = COS ( LAMDO-EET A ) 

SINDI = S IN ( LAMDO-EETA ) 

PSIDT = ( G * S I N ( F A ) * NLZ — VAR ( 8*3 ) * ( VWDZ*SlNDI-ETADZ*VW*COSDI 
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1 ) ) / ( V*C0S6A) 

TERM1 = VAR ( 8 * 3 ) * ( VWDZ*COSD I + VW*S I ND I *E TADZ ) 

VDOT = G*(NTZ-SINGA)+COSGA*TERMl 
GDOT = ( G* ( C OS ( F A ) *NL Z — COSG A ) —SINGA*TERM1 ) /V 

CL = <NLZ*WS)/Q 
Cl = (NTZ*WS)/Q 
63 G 1 = 57 .29578*GDOT 

MUDOT = <DADZ*SINGA) /COSMU - ( VDOT ) / ( A *COSMU*MSQ ) 

G2 = 57.29578*PSIDT 
G3 = 57.29578*MUDOT 

VG = SQRT(VAR( 8*1 ) *V AR ( 8*1 ) + V AR ( 8 * 2 ) *VAR ( 8 * 2 ) + V AR ( 8 * 3 ) *V AR ( 8 * 3 ) ) 

WRITE ( LUN I T * 4 ) VAR ( 8 * 1 ) * VAR ( 8 * 2 ) * VAR ( 8 > 3 ) *G1*G2»G3 

WRITE ( LUN I T * 5 ) Cl * CL *M * Q * V * VDOT 

WRITE ( LUN I T * 6 ) VWX»VWY»A»Vo 

XM = VAR (1*1) 

YM = VAR (1*2) 

ZINIT = VAR ( 1*3) 

TIME = TIND 
100 RETURN 

1 FORMAT ( 1 H 1 30X *MANEUVER DATA - OPTION TWO*) 

2 FORMAT ( *0T I ME =* Fll.l** H =*F11*1) 

3 FORMAT ( *0XG =* Fll.l** YG =*F11.1** ZG =*F11.1** GAMMA 

1 =* F 6 • 2 * * PSI =* F6 • 2 * * MU =*F6.2) 

4 FORMAT ( *OXGDOT =* Fll.l** YGDOT =*F11.1** ZGDOT =*F11.1** GAMDO 
IT =* F 6 • 2 * * PSIDOT =* F6.2** MUDOT =*F10.4) 

5 FORMAT ( *0CT-CD =* F11.6** CL =*F11.6** M =*F11.2** Q 

1=* F 8 • 1 * * V =*F8.1** VDOT =*F12.2) 

6 FORMAT ( *0UX =* F11.4** UY =*F11.4** A = *F 1 1 . 2 * 1 8 X * * VG 

1 =* F8.1) 

7 FORMAT ( 1 H 1 30X *MANEUVER DATA - OPTION ONE*) 

END 


SUBROUTINE INPOL(X>Y*T » FT ) 

C 

C GENERAL INTERPOLATION SUBROUTINE 
C 

C X IS THE ABSICISSA ARRAY 

C Y IS THE ORDINATE ARRAY 

C T IS THE ABSICISSA VALUE FOR WHICH A 

C CORRESPONDING ORDINATE VALUE * F T * WILL BE FOUND 

C 

c 

DIMENSION X( 2 5 ) * Y ( 2 5 ) 

IF(X( 1 ) — T) 11 *11*12 

11 DO 1 1=1*25 

IF < X( I )-T ) 1 *2 *3 

1 CONTINUE 

12 WR I T E ( 3*10) 

10 FORMAT ( 1X*29H ABSICISSA VALUE OUTSIDE RANGE) 

GO TO 4 

2 FT = Y ( I ) 

GO TO 4 

3 K= I — 1 

FT =Y(K) +<Y(I)-Y(K))*(T-X(K))/(X(I)-X(K)) 

4 RETURN 
END 
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SUBROUTINE LAGRA(T.N»TIME.TABLE»ANS1.ANS2»NSWIT) 
DIMENSION TIME! 10) , TABLE! 10) 

THIS ROUTINE FINDS THE INTERPOLATED VALUES OF A FUNCTION 
THE POINT T AND STORES THE ANSWER IN ANSI. ALSO THE 
AND STORES THE RESULT IN ANS2. 

NSW I T IS SET TO 0 FOR BOTH. 1 FOR ONLY FUNCTIONS -1 


NN = N— 1 
DO 3 L =2 *NN 
IF ( T- T I ME ( L+l ) 110.3.3 
3 CONTINUE 
L = N-l 

10 B1 = TIME(L-l) - TIME (L) 

B2 = TIME(L-l) - TIME (L+l) 

B3 = TIME(L)— TIME(L+1 ) 

A1 = T — T I ME ( L-l 1 
A2 = T - TIME(L) 

A3 = T - TIME! L+l > 

Cl = TABLE(L-1)/(B1*B2) 

C2 =-TABLE(L)/(Bl*B3) 

C3 = TABLE(L+1)/(B2*B3) 

I F ( NSW I T ) 18.15.15 

15 ANSI = A2*A3*C 1 + A1*A3*C2 + A1*A2*C3 
I F ( NSW I T ) 18.18.19 

18 ANS2 = C 1* ( A2+A3 ) + C2*(A1+A3) + C3*(A1 + A2 1 

19 RETURN 
END 


TABLE AT 
DERIVATIVE 

FOR ONLY DER 
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SUBROUTINE AREA 
REAL MU,NEW,MuDOT 

REAL M.LAMDO * MM ( 2 5 ) * LAMR < 2 5 ) » N T ( 2 5 ) . NL ( 2 5 ) » LA 
EXTERNAL OUT PU , DER I Z , PASS . DERPR ,OUTPR 
DIMENSION PHI (20) » T AU < 20 ) « F ( 2 0 » 20 ) , F I A ( 2 5 ) ♦ 

1 ALT ( 25 ) * TEMPI 25 ) » PRESS ( 2 5 ) *RHO< 2 5 ) »HW ( 25 > * 

2 VVW ( 2 5 ) * ETA < 2 5 ) » T T ( 2 5 ) »GR ( 2 5 ) 

DIMENSION C < 5 ) 

COMMON N*H*T » ACCUR , NODUM , REDUC, NODUB »ENDNO*ENDVL ,NSAVE » FLAG. Z» VAR 
1(14,8) 

COMMON IUNIT ,ISTAN,IALT ,1 TEMP , I RHO, I M * NPH I , NT AU , N ALTS * NW I ND , N I M 
1, WS, HG,H0, R E, PHI »TAU.F*FI A *ALT* TEMP .PRESS, RF, HW , V VW , E T A , MU , MUDOT , 
2TT »MM»GR»LAMR,NT *NL*LA*TINIT »GD0T ,DTPRN,DELZ1 , DELPT , F I MAX , DELF I ,PR 
3IN»MTIME*KTIME ,M,V ,VWX ,VWY *GAMR ,PSI DT ,RHOO,LAMDO ,TPRI N, T I Mt *Z INI T , 
5XM,YM,A,BETA, LUN IT ,DADZ*VW»EETA»VWDZ»ETADZ» MUNIT.FI , FA, DEL 
COMMON VWXX*VWYY»UN*UT*P*D»G0*G1»G4,G5 
COMMON SUM , I , DS I M , C , ZS I MP , I OUT , AGE , I PR I N , ZPR I N 
COMMON G2*G3»G50»G6*G7»G8»G9*G10*G11»G13»H1»H2»H3*H4»H5»TP»H7 
COMMON CF,SF.STH*CTH,THTAO* CU , SU , CG , SG , NE W , C N , SN , AO , UNO , 

1UT 0 , COO »CO»CT ,ST , I RH, ARA ,ML 
COMMON CONI ,CON2,CON3,CON4,L,PRES,DELP,DPMAX 
EQUIVALENCE ( PRESS! 1) , RHO ( 1 ) ) , ( Z , T I ND ) 

N0DUM=0 
REDUC = 0 . 

2000 CALL SSWT CH ( 6 , N ) 

GO TO ( 2001 ,2002 ) >N 

2001 CALL RESTA 

C INTEGRATION INITIALIZATION 

2002 KT I ME =K T I ME+ 1 

IF ( KT I ME- 400 ) 2003,2003,2004 

2004 TIME = T T < N I M ) 

GO TO 105 

2003 ML = 1 
N = 6 

NSAVE = 1 

ACCUR = 1. 

T = 1. 

NODUB = 1 
ENDNO=0. 

ENDVL=0. 

Z = Z I N I T 
VAR ( 1 , 1 ) =0. 

VAR ( 1 , 2 ) = 0 • 

VAR ( 1 ,3 ) =0. 

VAR (1,4) =XM 
VAR (1,5) =YM 
VAR ( 1 ,6 ) =T IME 

AREA AND RAY INITIALIZATION 

CF=CQS( FI*. 174532 9E-01 ) 

SF=SIN(FI*.1745329E-01 ) 

STH= ( BETA/M ) *CF*COS( GAMR ) -( 1 . /M ) *S I N ( GAMR ) 

CTH=SQRT ( 1 •- STH*STH ) 

ThTAO=ATAN(STH/SQRT( 1 • -ST H*ST H ) ) 

CU = COS (MU) 

SU = SIN (MU) 

CG = COS (GAMR) 

SG = SIN (GAMR) 

NEW = LAMD0— AT AN(CU*SF/(SU*CG+CU*CF*SG> ) 

CN *COS ( NEW ) 
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SN = S I N ( N E W ) 

TP IS A TEMPORARY LOACAT ION 

TO FIND THE SPEED OF SOUND AT INITIAL 2 
GO TO ( 1 *2 > * ISTAN 

1 CALL AT MO (Z*TP*P*D*AO) 

GO TO 3 

2 CALL I NPOL (ALT *TEMP*Z*TP) 

AO = 49 .020576*SQRT ( TP ) 

TO FIND WIND SPEED AT INITIAL Z 

3 S I NET = S I N ( EET A ) 

COSET=COS< EETA) 

C TO FIND DU/DZ 

TP =VWDZ*( SINET*SN + COSET*CN)-VW*ETADZ*( S I NET *CN-COSE T *SN ) 

UNO=VWX *SN+VWY *CN 

UT 0=—VWX *CN+VWY *SN 

COO= AO/CTH 

CO=CQO— UNO 

G2=1.+SG/(SU*STH) 

G3 =CG*CU*SF / CTH 

G50=RHOO*AO*AO 

CTH2=CTH*CTH 

G6=CU/CTH2 

G7 =CU*SG+SU*CG*CF 

G8 = STH/ (CTH*CTH2 ) 

G9=CU*STH*SF/CTH2 
G10=STH*SF 
G11=CG*SF/CTH2 
G1 3 =SU*CG+CU*SG*CF 

Hl=UTO*PSIDT-< G8*G7*A0-UT0*G1 1 ) *MUDOT- ( A0*G8*G1 3-UT0*G9 ) *GDOT 
1 -< AO*SG/SU)*(TP -DADZ/CTH) 

TP * CU*SF*CG 

H2=PSIDT+G1 1*MUDOT+G6*G10*GDOT 

H3=-UT0*G6*G7-A0*G8*TP 

H4=“G6*G7 

H5=G3*V+UT0*G2 

H7=Hl*H4-H2*H3 

CONI = CV*TP*G3)/ (STH*STH*AO) + (G6*G7*G2* CTH)/STH 
CONI = 1 •/( STH*CTH *SQRT < CON 1*G50*STH*CTH ) ) 

CON2 = (2.4*CONl>/CO 
SUM=0. 

I OUT = 1 
AGE=0. 

IPRIN = 1 
IRH * 2 
ZS I MP = Z I N I T 

C FIND ZPRIN THE FIRST PRINT OUT POINT 

DELP = ZINIT/DELPT 
IDEL = DELP 
TP = IDEL 

IF (TP-DELP) 558 *559 * 559 
559 IDEL = IDEL -1 
558 I F ( I DEL ) 551 *552*553 

551 IDEL = 0 

552 I OUT * 2 
ZPRIN = .01 
GO TO 183 

553 ZPR I N = I DEL*DELPT 

C FIND DSIM FOR INITIAL INTEGRATION STEP 
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183 DELP = ( Z I N I T - ZPRIN)/2. 

184 DELP =DELP/2. 

IF (DELP-DELZ1 ) 185*185,184 

185 DS IM=— DELP 
H = DS I M 

DELP = ( DS I M* . 4 ) /CO 

C< 1 ) = DE LP 
C( 2 ) = 4 • *DELP 
C( 3 > =2.*DELP 
C( 4 > =C ( 2 ) 

C< 5) =C( 1 ) 

1 = 0 

CALL NODEIDERIZ, PASS, OUTPU, PASS) 

GO TO (310,99,100) ,IRH 
99 GO TO ( 300 , 310 ) , IOUT 
300 IOUT = 2 
GO TO 184 

310 GO TO (30,40) ,ISTAN 
30 CALL ATMO( 0. , TP, PRES , D, A ) 

GO TO 50 

40 CALL I NPOL (ALT , PRESS ,0. , PRES ) 

C PRESSURE CALCULATIONS 
50 N = 1 
T=3. 

ACCU R = 0 . 

NODUB = 0 
NSAVE = 1 

TP = ABS ( FI*, 1745329E— 01-FA ) 

4 IF ( PHI ( ML ) -TP ) 6,5,5 

6 IF (ML-NPHI) 7,5,5 

7 IF ( ABS ( PH I (ML)-TP) —ABS (PHI (ML+ll-TP) ) 5,5,8 

8 ML = ML + 1 

GO TO 4 

5 CONI = SORT ( ( G50*CU*CTH2 ) /SU ) 

CON2 = SU/(CO*CTH) 

CON3 = • 5*AGE 

C0N4 = SORT ( G5/ ( G4*ARA ) ) 

TIND = T AU ( 1 ) *LA 
DPMAX = 0. 

VAR ( 1,1) = 0. 

VAR (8,1) = 0. 

VAR (1,2) = F ( 1 »ML ) 

ENDNO = 0. 

ENDVL = TIND 
L = 1 

CALL OUTPR 
DO 11 L = 2 , NT AU 
FLAG = 0. 

ENDVL = T AU ( L ) *L A 
ENDNO = 4. 

CALL NODE ( DERPR , DERPR , OUTPR, PASS) 

11 CONTINUE 

GO TO ( 110,400,110) , I RH 
110 KTIME = KTIME-1 
GO TO 100 

400 WRITE (MUNIT) F I , VAR ( 1 , 4 ) , VAR ( 1 , 5 i , VAR ( 1 , 6 ) , DPMAX 

100 F I =F I +DEL 

IF ( ABS ( F I ) -F IMAX ) 2000,2000,101 

101 IF(FI) 105,102,102 

102 FI = -DELFI 
DEL = -DEL 
GO TO 2000 

105 RETURN 
END 
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SUBROUTINE DERIZ 
REAL MU » NEW * MUDOT 

REAL M.LAMDO ,MM ( 2 5 ) , LAMR ( 2 5 > . NT ( 25 ) » NL ( 2 5 ) . LA 
DIMENSION PHI (20) , TAU ( 20 ) *F(20»20) ♦FIA125) * 

1 ALT (25) .TEMP (25) .PRESS (25) .RHO<25) .HW(25) . 

2 VVWI 25) *ETA(25) . T T ( 2 5 ) . GR ( 2 5 ) 

DIMENSION C ( 5 ) 

COMMON N. H. T. ACCUR . NO DUM. R EDUC. NODUB. EN DNO. END VL . NS AVE . FLAG. Z. VAR 
1(14.8) 

COMMON I UN IT. I STAN. I ALT * I TEMP . I RHO. IM.NPHI . NT Au » NALTS . NW I ND . N I M 
1 .WS.HG.HO.RE.PHI . T AU . F . F I A . AL T . TEMP . PRESS . RF . HW ,VVW . ETA .MU . MUDOT . 
2TT .MM.GR . L AMR .NT .NL.LA.T INIT »GDOT .DTPRN.DEL21 .DELPT .F I MAX .DELE I .PR 
3IN.MTIME.KTIME.M.V.VWX.VWY.GAMR.PSIDT . RHOO .LAMDO.TPRIN.TIME.ZINIT. 
5XM.YM.A.BETA. LUN IT .DADZ.VW.EETA.VWDZ.ETADZ. MUNIT.FI .FA. DEL 
COMMON VWXX.VWYY.UN.UT .P.D.G0.G1 .G4.G5 
COMMON SUM .1 * DS I M . C . ZS I MP . I OUT . AGE . I PR I N . 2PR I N 
COMMON G2.G3.G50.G6.G7.G8.G9.G10.G11.G13.H1.H2.H3.H4.H5.TP.H7 
COMMON CF.SF.STH.CTH.THTAO. CU » SU »CG . SG . NEW . CN . SN » AO . UNO . 

1UT0 .COO.CO.CT.ST. IRH.ARA.ML 
EQUIVALENCE ( PRESS ( 1 ) ,RHO ( 1 ) ) , ( AST . VV ) . ( SC . TP ) 

CALL SSWTCH ( 5 . J ) 

GO TO ( 100.4) . J 

100 WRITE ( LUN IT .155) Z . V AR ( 1 . 4 ) ♦ VAR ( 1 . 5 ) . VAR ( 1 . 6 ) . V AR ( 1 . 1 ) . VAR ( 1 . 2 ) . V 
1ARI1.3I 

155 FORMAT ( 3H Z=.E12.4,3H X=.E12.4.3H Y=.E12.4.4H T= E12.4.4H I1=.E12 
1.4.4H 12= E14.4.4H 13= .E14.4) 

C CALCULATE SPEED OF SOUND . DENS I TY . PRESSURE . AND WIND SPEED 

4 I F ( Z ) 5.6.6 

5 I F ( Z-H* .01) 10.10.7 

7 Z = 0. 

6 GO T 0 ( 1,2) .ISTAN 

C TP IS A TEMPORARY LOCATION 

1 CALL ATMO(Z.TP.P.D.A) 

GO TO 3 

2 CALL INPOLI ALT . TEMP » Z . TP ) 

A=49. 020576*SQRT ( TP ) 

CALL INPOLt ALT, PRESS. Z.P) 

GO T0( 11 » 12 ) , IRHO 

11 D=P/ ( 1716. *TP ) 

GO TO 3 

12 D = P 

P=D*< 1716. *TP) 

3 CALL LAGRA(Z.NWIND»HW.ETA»TP»UT .1 ) 

CALL LAGRA(Z»NWIND.HW, VVW .VV.UN.l ) 

VWX X = VV*S I N ( TP) 

VWYY=VV*COS ( TP ) 

UN=VWXX*SN+VWYY*CN 
UT=-VWXX*CN+VWYY*SN 
CT=A/ (CO+UN ) 

ST = SQRT ( l.-CT*CT ) 

AST=A*ST 

SC=ST/CT 

GO=CT/ST 

G1=G0*G0*G0/ ( A*A ) 


VAR (8*1) 

= -Gl 

VAR (8*2) 

=-Gl*UT 

VAR ( 8*3) 

=-Gl*UT*UT-GO 

VAR (8*4) 

= -SN/SC+VWXX/AST 

VAR (8*5) 

=— CN/SC+VWYY/ AST 

VAR (8*6) 

=— 1 • /AST 

VAR ( 1 *7 

) = CT 

VAR ( 1*8 

) s D*A*A 


C 

10 RETURN 
END 
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SUBROUTINE OUTPU 
REAL MU * NEW ♦ MUDOT 

REAL M.LAMDO » MM ( 2 5 ) » L AMR ( 2 5 ) » NT ( 2 5 ) , NL ( 2 5 ) » L A 
DIMENSION PHI ( 20 ) »TAU< 20) *F ( 20 » 20 ) »FI A ( 25 ) ♦ 

1 ALT(25).TEMP(25).PRESS(25).HW(25) »VVW (25)*ETA(25)*TT(25)*GR(25) 
DIMENSION C < 5 ) 

COMMON N »H* T » ACCUR .NODUM .REDUC . NODUB » ENDNO » ENDVL .NSAVE » FLAG* Z » VAR 
1(14.8) 

COMMON I UN IT. I STAN. I ALT. I TEMP, I RHO. I M . NPH I » NT AU . NALTS » NW I ND » N I M 
1 .WS.HG.HO.RE.PHI » T AU . F » F I A . ALT , TEMP , PRESS . RF . HW .VVW. » ETA » MU * MUDOT . 
2TT.MM,GR,LAMR.NT»NL.LA.TINIT.GD0T»DTPRN,DELZ1 .DELPT .FIMAX.DELFI .PR 
3IN.MTIME.KTIME.M.V , VWX . VWY » GAMR » PS I DT » RHOO » LAMDO » T PR I N .TIME.ZINIT. 
5XM,YM,A.BETA,LUNIT,DADZ.VW.EETA.VWDZ,ETADZ. MUNIT »FI » FA. DEL 
COMMON VWXX . VWYY »UN »UT .P.D.GO.Gl .G4 .65 
COMMON SUM.I.DSIM. C , ZS I MP . I OUT ♦ AGE , I PR I N » ZPR I N 

COMMON G2,G3»G50»G6»G7,G8.G9,G10,G11,G13.H1,H2,H3»H4.H5.TP»H7 
COMMON CF .SF »STH .CTH » THTAO , CU . SU »CG » SG » NEW . CN , SN » AO » UNO . 

1UTO.COO.CO.CT.ST. IRH.ARA.ML 
COMMON CONI » C0N2 »C0N3 » C0N4 » L » PRES » DELP 
CT = VAR (1,7) 

IF.(VAR( 1,7)-. 9994) 200.210,210 

210 FLAG = -1. 

IRH = 1 

GO TO 6 

200 IF (Z/ZSIMP-1. 00001 ) 1 ,1.300 

1 1=1+1 

C CALCULATE AREA 

ST = SORT ( 1 .-CT*CT ) 

G4=ST*CT 
G5 = VAR (1.8) 

ARA=H3*( VAR ( 1 , 2 ) *G2-VAR (1,1 ) *H5 ) -H4* ( VAR (1.3) *G2-VAR (1.2) *H5 ) 

1 +H7* ( VAR ( 1 » 1 ) *VAR ( 1 ,3 )-VAR (1»2)*VAR(1,2>) 

C CALCULATE INTEGRAND OF AGE FUNCTION AND SUM 
IF (ARA) 410.410.415 
410 TAUA = 0. 

IF (Z-ZINIT) 420.350,350 
350 TP = NEW/. 1745329E-01 

WRITEfLUNIT .556) FI, TP .COO.CO 
556 FORMAT ( //30X ,*RAY .AREA AND PRESSURE CALCULATIONS*// 

110X»3HFI=,F10.3.10X»3HNU=»F10.3»10X»4HC00=»F10.3,7X»3HC0=»F10.3// 
26X.1HX.15X.1HY.15X.1HZ.15X.1HT.12X.4HAREA.11X. 3HAGE .10X.9HC0S THET 
1A) 

GO TO 50 
430 IRH = 4 

GO TO 421 

420 IRH = 3 

421 FLAG = -1. 

AGE = — SUM+C0N2*SQRT ( Z I N I T— Z ) 

GO TO 6 

415 TP = SORT ( Z I N I T— Z ) 

I F ( G4*G 5 ) 430,430,427 

427 TAUA= 1 * /SORT ( G5*ARA*G4*G4*G4 ) -CON 1 / TP 

416 SUM= SUM+Cl I )*TAUA 
401 IF(I-5)50.5»50 

5 AGE = -SUM+C0N2*TP 
SUM=C ( 1 )*TAUA+ SUM 
1 = 1 

HXS=VAR (1,4) 

HYS = V AR (1,5) 

HZ S = Z 

HT S*VAR (1.6) 
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AREAS=ARA 

CTS=CT 

C CHECK TO SEE IF AT ZERO OR P&IN OR PRINT INTERVAL 
IF (Z/ZPP I ^-1.00001 ) 66*66*9 
21 IPRIN = 2 
6 V A R ( 1 • 4 ) = H X S 
VAR ( 1 * 5 ) -H YS 
Z = HZS 

VAR( 1 *6)=HTS 
ARA= AREAS 
CT=CTS 

C SAVE VARIABLES FOR USE IF THETA OR AREA APPROACH ZERO 

C RESTORE VARIABLES WHEN THETA OR AREA APPROACH ZERO 

C 

c 

WRITE ( LUNIT * 61 ) VAR (1*4) * VAR ( 1 * 5 ) * Z * VAR ( 1 * 6 ) * ARA * AGE *CT 
61 F0RMAT(1X*7(E13*5*?X) ) 

68 IF(Z-.l) 10*10*100 
10 FLAG = -1. 

GO TO (88*300*77*99) * IRH 
99 WRITE(LUNIT*90) 

90 FORMAT ( /20X*27HNEGATIVE SQUARE ROOT IN AGE ) 

IRH = 3 
RETURN 

77 WRITE(LUNIT *70) 

70 FORMAT ( /20X >9HAREA ZERO) 

RETURN 

88 WRITE (LUNIT *80) 

80 FORMAT ( /20X * 14HRAY HORIZONTAL) 

RETURN 

9 GO TO ( 8 ♦ 50 ) * IPR IN 
8 IF(Z-PRIN) 21*21*50 

66 ZPRIN = ZPRIN -DELPT 

IF ( ZPRIN ) 64*64*65 

64 ZPRIN = .1 

65 I F ( Z— PRIN) 67*67*6 

67 IPRIN = 2 
GO TO 6 

100 GO TO ( 14*50 ) * I OUT 
14 DELP=DFLPT / 2 • 

GO TO 10 

50 ZSIMP=7SI M P+DSIM 

GO TO (88 *52 *77) » IRH 
52 IF(ZSIMP-1.) 51*51*300 

C ADJUST ZSIMP SO THAT DIVIDE BY ZERO DOES NOT OCCUR 

51 ZSIMP = .1 
300 RFTURN 

END 
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SUBROUTINE DERPR 

REAL M.LAMDO »MM ( 2 5 ) * LAMR (25) »NT(25)»NL(25) * LA » MU .MUDOT *NEW 
DIMENSION PHI (20) »TAU<20) .F(20.20).FIA(25>»C(5>, 

1 ALT ( 25 ) .TEMP (25).PRESS(25).HW(25) »VVW (25).ETA(25).TT(25)*GR(25) 
COMMON NNODE.HNODE.TNODE.ACC.NODUM.BYHAL . NODUB. ENDNO » ENDVL » NSAVE » 
5FLAG.TIND.VAR( 14.8 ) 

COMMON IUNIT. I STAN . I ALT . I T EMP . IRHO. I M . NPH I » NT AU . NALTS » NW I ND . N I M 
1. WS.HG. HO, RE. PHI .TAU.F.FIA.ALT .TEMP. PRESS. RF. HW , VVW . ETA . MU .MUDOT » 
2TT.MM.GR. LAMR .NT.NL.LA.TINIT »GDOT .DTPRN.DELZl , DELPT . F I MAX , DELF I » PR 
3 IN .MTIME.KT I ME .M.V.VWX.VWY ,G AMR. PS I DT , RHOO .LAMDO » TPR I N. TIME.2INI T . 
5XM.YM.A.BETA. LUN I T .DADZ.VW.EETA . VWDZ , ETADZ » MUN I T ,FI »FA .DEL 
COMMON VWXX.VWYY .UN.UT .P.D.G0.G1. G4 . G5 
COMMON SUM.I.DSIM. C , ZS I MP » I OUT » AGE » I PR I N » ZPR I N 

COMMON G2.G3.G50.G6.G7.G8.G9.G10.G11.G13.H1.H2.H3.H4.H5.TP.H7 
COMMON CF.SF.STH.CTH. THTAO.CU.SU, CG.SG.NEW.CN.SN.AO.UNO.UTO. COO . CO 
l.CT.ST.IRH.ARA.ML 
COMMON CONI .CON2.CON3 »C0N4 » L 
TP = TIND/LA 

CALL INPOLf TAU.F ( 1 .ML ) .TP.VARf 1.2)) 

VAR (8.1) = C0N1*VAR( 1 ,2 ) 

CALL SSWTCH ( 5 » J ) 

GO TO ( 100.1 ) » J 

100 WRITEfLUNIT.155) T I ND . VAR ( 1 , 1 ) , VAR ( 8 , 1 ) 

155 FORMAT ( 3H Z= E 1 5 .4 , 9H I NT EGR AL= E 14 . 5 » 1 OH I NT EGR AND= E15.4) 

1 RETURN 
END 

SUBROUTINE OUTPR 

REAL M.LAMDO . MM ( 2 5 ) . LAMR ( 2 5 ) » NT ( 2 5 ) . NL ( 2 5 ) » LA , MU » MUDOT .NEW 
DIMENSION PHI (20) .TAU(20) ,F(20.20)»FIA(25).C(5), 

1 ALT (25) »TEMP(25) »PRESS(25) , HW ( 25 ) , VVW (25>»ETA(25) » T T < 25) »GR(25) 
COMMON NNODE .HNODE .T NODE. ACC, NODUM.BYHAL .NODUB » ENDNO , ENDVL . NSAVE » 
5FLAG,TIND»VAR( 14,8) 

COMMON IUNIT .ISTAN.IALT. I TEMP. IRHO. I M, NPH I . NT AU .NALTS ♦ NW I ND » N I M 
1. WS.HG, HO, RE .PHI . T AU . F , F I A » AL T » TEMP , PRESS , RF » HW , VVW . ET A , MU . MUDOT » 
2TT, MM, GR .LAMR. NT .NL.LA.T IN IT. GDOT » DTPRN.DELZl .DELPT » F I MAX , DELF I .PR 
3 IN.MT I ME ,KT I ME »M.V .VWX »VWY »GAMR ,PSI DT » RHOO » LAMDO , TPR IN. TIME. ZI NIT. 
5XM » YM » A .BETA » LUN IT.DADZ.VW.EETA.VWDZ.ETADZ, MUN IT.FI.FA. DEL 
COMMON VWXX.VWYY .UN.UT » P , D .GO » 61 . G4 , G5 
COMMON SUM.I.DSIM, C » ZS I MP , I OUT , AGE » I PR I N . ZPR I N 

COMMON G2»G3»G50»G6»G7.G8»G9»G10.G11,G13,H1»H2,H3»H4»H5»TP»H7 
COMMON CF.SF .STH.CTH. THT AO »CU .SU.CG.SG.NEW.CN »SN » AO .UNO »UT0» COO. CO 
1 »CT .ST » IRH.ARA.ML 

COMMON CONI » C0N2 ,C0N3 »C0N4 , L , PRES »DELP .DPMAX 
IF(L-l) 5,5.6 

5 WRITE ( LUN IT .200) 

200 FORMAT ( //5X ♦ 1HL » 10X » 1HF ♦ 1 IX, 5HVE( L ) »7X »4HX I 0.8X.4HXI 1.8X.5HS INT 

1.7X.1HS.11X.2HP1.10X.5HDEL P> 

6 IF ( ENDNO- 4. ) 10.20,20 

10 XIO = TIND*C0N2 

XII = X 1 0— AGE* VAR (8.1) 

SINT = C0N2*VAR< 1.1) 

S = SINT -CON3*VAR(8,l)*VAR(8»l> 

DELP = VAR( 8.1 )*C0N4*RF 
IF(ABS(DELP)-DPMAX) 50.50,60 
60 DPMAX = ABS(DELP) 

50 PI = DELP/PRES 

WRITE ( LUNIT .2 ) TIND.VAR < 1 .2 ) » VAR( 8 » 1 ) .XIO.XIl . S I NT » S , P 1 » DELP 

2 FORMAT (9F12.6) 

IF ( ENDNO) 30.30.20 
30 FLAG = -1. 

20 RETURN 
END 
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SUBROUTINE NODE(COMPP »CQMP Y » COMPT » COM PE ) 

DIMENSION VAR ( 14*1 ) * A ( 4 ) * B ( 4 ) * C ( 4 ) 

COMMON N tH*T *PHSIG .NODUMt BYHAL * NODUB * ENDNO * END VA * NSAVE * FLAG » X * 
1 VAR 

C INITIALIZE 

MODF ( XSAVE * PERR ) = XSAVE -FLOAT ( I F I X ( XSAVE/PERR ) ) *PERR 
FLAG= 0 • 

500 IF (BYHAL >502*501*502 

501 BYHAL = .& 

502 IF ( ENDNO) 503*504*503 

C ENDPOINT COMPUTES H 

503 H= ( ENDVA -X)/ENDNO 

C PREPARE FOR RKG 

504 CALL COMPD 

505 CALL COMPT 

IF ( FLAG) 560*506*506 

506 XSAVE=X 

C PROGRAM USING RKG AS A STARTING PROCEDURE 

I F ( A ( 1 ) — • 5 ) 400 >40 1 * 400 

400 A ( 1 ) = • 5 

A< 2) = . 29289322 
A ( 3 ) = 1 • 7071 068 
A ( 4) =.16666667 
B ( 1 ) = 1 • 

B ( 2)=. 29289322 
B ( 3)=1. 7071068 
B(4)=. 33333333 
C( 1 ) = • 5 

C( 2) =.29289322 
C( 3)=1. 7071068 
C ( 4 i = • 5 

401 DO 402 I = 1 * N 

402 VAR (6*1 ) =0 • 

J= 4 

GO TO 410 

403 DO 407 K = 1 * 4 
DO 404 I = 1 * N 
CK=H*VAR (8*1 ) 

R= ( A ( K ) *CK ) - ( B ( K ) *VAR (6*1)) 

VAR ( 1 t I ) = VAR ( 1 ♦ I )+R 

404 VAR (6*1 ) = VAR (6*1 ) +( 3 • *R ) - ( C ( K ) *CK ) 

IF ( K— 1 )405*405*413 

413 I F ( K-3 ) 406*405*406 
C NEW VALUE OF X 

405 X=X+ ( H/2 • ) 

CALL COMPD 
GO TO 407 

406 CALL COMPY 

407 CONTINUE 

IF ( NODUM )410*412*411 
412 NODUMa-1 

410 NN=N+NSAVE 

DO 408 1=1 *NN 
VAR(J + 1*I )=VAR( 1*1 ) 

408 VAR( J + 8*I )=VAR(8*I ) 

J= J-l 

IF (J) 409*409*403 

411 CALL COMPT 

IF (FLAG)560*410*410 

409 NSWHF= 1 

IF ( ENDNO) 507 *508 *507 
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c 


c 


c 


507 ENDNO=ENDNO-3. 

508 M= 3 

509 FLAG=0. 

510 X = X + H 


PROGRAM FOR THE PREDICTOR 

DO 450 I = 1 * N . , „ 

450 VAR< 1 ♦ I ) = ( 1 . 54765 11*VAR< 2. I ) )-( 1 • 8675052*VAR (3*1 ) ) + (2 
1 VAR (4. 1 ) )-( . 697352 8*VAR< 5.1) )+H*( ( 2 . 0022473*VAR < 9 ♦ I ) )■ 
2 VAR ( 10. I ) ) + ( 1 .8186108*VAR< 11. I) )-(.71432005*VAR(12*I> 
512 CALL COMPD 
PERR=0. 

PROGRAM FOR THE CORRECTOR 
DO 46 2 I = 1 . N 

TEMP = VAR(2»I ) +H* ( I .375*VAR( 8. I ) ) + ( • 79166667 *VAR (9. 1 ) ) 
l-( .2083333 3 *VAR ( 10 ♦ I ) )+( .041666667*VAR( 11. I ) ) > 


513 


460 


IF (PHSIG-1.) 463.464.463 
463 TEMPA=ABS ( ( TEMP-VAR ( 1 ♦ I ) > / T EMP ) 
GO TO 465 


464 TEMPA=ABS (TEMP-VAR (l.I)) 

465 VAR ( 1*1) =TEMP 

IF ( PERR-TEMPA )461 .462 .462 

461 PERR=TEMPA 

462 CONTINUE 

515 CALL COMPY 

516 IF < PERR-16. 219659*10. **(-T )> 517.517 .535 
NO HALVING NECESSARY 


517 NSWHF=0 

IF ( NODUM 1550.518.518 

518 IF (ENDN01519.520.519 

519 ENDNO=ENDNO-l. 

520 CALL COMPT 

IF ( FLAG ) 560.521 .521 
C IS DOUBLING POSSIBLE 

521 IF (PERR-(( 16.219659* 10 . ** ( -T ) ) /200 • ) ) 525.525.522 

522 M= 3 
528 J= 1 3 

523 DO 524 1*1. N 

524 VAR( J + l ♦ I ) =VAR ( J. I ) 


J= J-l 

IF (J) 509.509.523 

C DOUBLING 

525 M=M— 1 

526 IF < M ) 530.527.528 

527 IF (ENDNO)530»531*530 

530 IF (MODFIENDNO.2. ) 1528.531.528 

531 FLAG = 2 . 

CALL COMPE 

IF (FLAGI560.532.532 

532 IF(NODUB)522.529.522 
529 NN=N+NSAVE 

DO 533 1=1. NN 
VAR (2.1 ) =VAR < 1 . 1 ) 

VAR (4.1 ) =VAR (5.1 ) 

VAR (5.1 ) =VAR (7.1 ) 

VAR ( 9.1 ) = VAR <8.1 ) 

VAR ( 11 . I ) =VAR ( 12.1 ) 

533 VAR ( 12.1 >=VAR< 14.1 ) 

H=2.*H 

IF (ENDN01534.508.534 

534 ENDNO=ENDNO/2 . 

GO TO 508 


.0172069* 
-12.0316877* 
) ) 
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C HALVING 

535 FLAG=ABS ( BYHAL ) 

CALL COMPE 

565 IF (FLAG)560*561»561 
561 IF ( NODUM 1537.537*536 

536 CALL COMPT 

IF (FLAG)560*537»537 

537 IF (BYHAL -1 . ) 548 . 5 1 7 . 5 1 7 

548 I F ( ENDNO - 1 • ) 542*563*543 
563 ENDNO= 2.*ENDNO 

H=H/2. 

FLAG= . 5 
CALL COMPE 

IF (FLAG) 560*543*543 

543 ENDNO=ENDNO/ BYHAL 
NN=N+NSAVE 

542 IF <NSWHF)538.540.538 
C REPEATED HALVING 

538 DO 539 I=1*NN 
VAR ( 1 » I ) =VAR (5*1 ) 

539 VAR (8*1 )=VAR( 12.1 ) 

X=XS AVE 

IF (ENDNO) 549.549 *544 

544 ENDNO=ENDNO+3. /BYHAL 

549 H=H*ABS (BYHAL ) 

GO TO 506 

540 DO 541 1*1. NN 
VAR (1*1 ) =VAR (2.1 ) 

541 VAR (8*1 ) =VAR (9.1 ) 

X = X-H 

GO TO 549 

C DUMMY OUTPUTTING 

550 X=XS AVE+H 

IF (ENDNO)551.552.551 

551 ENDNO=ENDNO+2. 

552 K*3 
NN=N+NSAVE 

DO 553 1 = 1 * NN 
VAR (6.1 ) =VAR ( 1 . I ) 

553 VAR(13*I)=VAR(8.I ) 

557 DO 554 1=1. NN 

VAR ( l.I )=VAR(K+1»I ) 

554 VAR (8*1) =VAR ( K + 8 . I ) 

CALL COMPT 

IF (FLAG)560*562.562 
560 RETURN 
562 X=X+H 
K=K-1 

IF ( K ) 568.558.555 

555 I F ( ENDNO >556.557.556 

556 ENDNO*ENDNO-l. 

GO TO 557 

558 DO 559 1=1. NN 
VAR ( 1 . I ) =VAR (6.1 ) 

559 VAR (8.1 ) =VAR (13.1 ) 

N0DUM*0 

GO TO 518 
END 
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non 


SUBROUTINE PASS 

RETURN 

END 


C 

C 


SUBROUTINE 
ROUTINE FOR 


G ROUN 

INTERPOLATION 


OE GROUND INTERSECTIONS 


NRMAX.HN.TN.AC.MUNI T.BY.LUNI T.T (400) .XI 400) . Y (400) .P1400) 


OF SAVED GROUND INTERSECTIONS) 


DISK OR TAPE 
FI.X(I).Y(t).T(I).P(I> 


COMMON 
1N( 41 ) 

I = 1 
K = 1 
NPHSW = 1 
WRITE (LUNIT.667) 

667 FORMAT ( 1H1 .20X .34HLIST 
DO 150 NREC = l.NRMAX 
READING A RECORD FROM THE 
30 READ ( MUN I T ) 

WRITE (LUNIT.666) F I . X ( I ) . Y ( I ) ♦ T ( I ) »P ( I ) » I 
666 FORMAT ( 9F20.5 .2110) 

IF(FI) 50.200.100 

50 GO TO (51.100) .NMNSW 
FIRST NEGATIVE ANGLE 

51 NMNSW = 2 

CALL MOVE! I . 1+1) 

CALL MOVEt IMIN.I ) 

N ( K ) = I 

K = K + 1 
1 = 1+1 

FIND THE TIME MINIMUM FOR CURVE 
TMIN ) 110.150.150 


100 

I F ( T ( 

I ) - 

110 

TMIN : 

= T ( I 


J = I 



GO TO 

150 

F 

I IS 1 

EQUAL 

2 00 

GO TO 

(201 

201 

NPHSW 

= 2 

220 

N( K) 

= I 


K = K 

+ 1 


IMIN 

= I 


NMNSW 

= 1 


GO TO 

110 


TO ZERO 
.210) .NPHSW 


CHECK IF TMIN IS AT FI = 0 
210 GO TO ( 211 .215 ) .NMNSW 
K I = 1 + 1 
CALL MOVEt I *KI ) 

CALL MOVE! IMIN.I ) 

N( K) = I 
K = K + 1 
I = 1+1 

IF(J-IMIN) 1000*220.230 
ISAVE = 1+1 
JMAX = N(K-1)-1 
LOOPS = IMIN +1 


211 


215 

230 


MOVE TABLES OUT OF PERMUTATION AREA 
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DO 231 KI = LOOPS. JMAX 
CALL MOVE ( K I * I SAVE ) 

231 ISAVE = ISAVE + 1 
C MOVING THE TABLES AROUND 

IF ( J-N ( K— 1 ) ) 240 .1000.250 
240 KK = J-IMIN-1 

LOOPS * I + UK +1 
IC = 1 

N< K-l ) = JMAX - KK 
DO 244 KI = IMIN.JMAX 
CALL MOVE ( LOOPS .K I ) 

LOOPS = LOOPS + IC 
I F ( LOOPS- ISAVE) 244.242.242 
242 LOOPS = I + KK +1 
IC = -1 
244 CONTINUE 
GO TO 220 
250 JMAX = J-l 
IC = -1 
LOOPS = J 

DO 254 K I = I M I N .JMAX 
CALL MOVE ( LOOPS .K I ) 

LOOPS = LOOPS + IC 
I F ( LOOPS - NtK-11) 252*254.254 

252 LOOPS = I +1 
IC = 1 

254 CONTINUE 

N( K-l ) = J 

GO TO 220 
150 I » I + 1 
KMAX » K-2 
WRITE! LUNIT .413) 

DO 450 K » 3. KMAX. 2 
WRITE (LUNIT. 412) 

KN = N ( K ) 

TVAL = T ( KN ) 

KI = K + 1 
KK = 1 
IC = 2 

400 CALL INTEP ( TVAL. XVAL .YVAL .PVAL.KK .KERR ) 

GO TO (405.499) .KERR 

405 WR ITE ( LUNIT .410 ) T VAL ♦ XVAL . YVAL .PVAL 

410 FORMAT ( 4F20 • 5 ) 

406 KK = KK +IC 
IF(KK-KI) 420.450.411 

420 IF(KK) 450.450.400 

411 KK = K- 1 
IC = -2 

GO TO 400 

222 FORMAT ( * NO DATA POINT ON CURVE * 15) 

499 WRITEtLUNIT.222) KK 
GO TO 406 
450 CONTINUE 
20 WRITE (LUNIT. 60) 

60 FORMAT ( //20X.1 7HPROBLEM COMPLETED) 

RETURN 

412 FORMAT ( / 10X * 4HT I ME ♦ 1 6X . 1HX . 1 9X . 1H Y . 19X . 8HPRESSURE ) 

413 FORMAT ( 1H1 .20X .26HGROUND INTERSECTION CURVES) 

1000 WRITEtLUNIT .1001 ) 

1001 FORMAT ( 13H KULSRUD GOOF) 

CALL EXIT 

END 


140 



SUBROUTINE I NTEP ( T VAL ♦ X ANS * YANS * PANS * KK * NERR ) 

C AN INTERPOLATION ROUTINE IN THE SPECIAL TABLES 

COMMON NRMAX*HN»TN*AC*MUNIT *BY*LUNIT ♦ T (400) *X(400) *Y(400) tP( 400) * 
1N( 41 ) 

IN = NOCK) 

INF = N(KK+1)-1 
DO 200 IS = INtINF 
I F ( T ( ISJ-TVAL) 200 * 100 » 300 
100 PERCT = 1. 

GO TO 302 

200 CONTINUE 

201 NERR = 2 
RETURN 

300 IFIIS-IN) 2 01* 201 *301 

301 PERCT = (TVAL-T(IS-l) )/CT( IS)-T<IS-1) ) 

302 XANS = PERCT*(X( IS)-X( IS-1 ) )+ X(IS-l) 

YANS = PERCT*(Y< IS)-Y( IS-l) >+ Y(IS-l) 

PANS = PERCT*(P( IS)-P( IS-1 ) ) + PCIS-l) 

NERR = 1 

RETURN 

END 


SUBROUTINE MOVE(KO) 

COMMON NRMAX*HN*TN#AC*MUNIT*BY*LUNIT*T<400> *X(400> *Y<400> *P(400) * 


IN ( 41 ) 
T( J) = 

Tl 

[<) 

X( J) = 

X 1 

IK) 

Y( J) * 

Y! 

IK) 

P( J) = 

PI 

(K) 

RETURN 

END 




SUBROUTINE RESTA 
COMMON C(1000) 

REWIND 9 
WRITE (9) C 
REWIND 9 
WRITE (6*1) 

1 FORMAT ( /23H0TERMINATED BY OPERATOR 
CALL EXIT 

2 RETURN 
END 


SUBROUTINE INITAL 

COMMON C(1000) 

REWIND 9 

READ (9) C 

REWIND 9 

RETURN 

END 
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MANEUVER DATA - OPTION ONE 
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